MATLAB: MEX matrix division gives different result

Question

I've used MATLAB's coder tool to create a MEX version of the matrix exponential function, to be used in another set of functions. The issue is, the MEX version gives different results than the original m-file.

After debugging it, I believe that the reason this is, is because the MEX file and the m-file do not do matrix division (\) the same. Or the MEX file has issues with it in the first place. All the variables leading up to the line where the matrix division occurs are equivalent on both sides.

This is the line where the issue occurs:

F = (V-U)\(2*U) + I

Where I is the identity matrix of the size of V and U.

What is the reason behind the discrepancy when a MEX file does matrix division, and how can I fix this issue? Can this line of code be re-written without the division?

Answer 1

I have no problem generating C code from such an operation.

Here is a test function I tried:

myfcn.m

function F = myfcn(U,V)
    I = eye(size(U));
    F = (V-U)\(2*U) + I;
end

here's a test script we'll use to verify the results:

test_myfcn.m

U = rand(5);
V = rand(5);
F = myfcn(U,V);

I start by launching the code generation tool (ccoder), create a new project set to produce a MEX-file, then add the myfnc.m function from before as entry-point. Then I define both input variables types as:

double (:Inf x :Inf)

which specifies an MxN matrix of unbounded size of type double.

Finally we can build the project. This produces myfcn_mex.mexw64.

Testing both the original M-function and the generated MEX-function, I get pretty much the same result (difference is in the order of machine epsilon):

>> F = myfnc(U,V);
>> FF = myfcn_mex(U,V);
>> norm(F-FF)
ans =
   1.4391e-14

Answer 2

MATLAB made a minor change to the EXPM algorithm in the 2014a release. MATLAB Coder implementations are separate, and the corresponding change has not yet been made in the algorithm for code generation, so some discrepancies are possible there. The new implementation does (V - U)\(2*U) + I, whereas the old one does (V - U)\(V + U). These are mathematically equivalent but will give different round-off behaviors in general.

AFAIK, there is no systematic difference in the quality of solutions to linear systems in MATLAB versus MATLAB Coder. The core algorithms are essentially equivalent, with round-off differences expected to creep in from various obscure sources. The residual may be smaller for either MATLAB or MATLAB Coder in a given case. If the difference in the solutions is large, it indicates that the problem being solved is ill-conditioned. I can explain more about this if you like, but it is covered in every numerical analysis textbook. Can you provide a concrete example? At least can you find out what cond(V - U) returns there when your problem is solved in MATLAB?