MEX usage of mldivide() produces incorrect/different results

I have used the coder toolkit to create a MEX'ed version of the matrix exponential function. Everything up until the point of the matrix division is equivalent. However, once it arrives at this line of code:

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

The result of the function becomes inaccurate. Although the values are small, the percent change from the expected (non-mex) answer can reach up to over 100%.

Why is this discrepancy happening, and what can I do to fix this so that, while the answers don't have to be equivalent, the percent change from the actual answer must be at an acceptable level?

Or perhaps there is another method to solve the equation without the use of "\". Note: I have already tried the inv(A)*B method.

*I have provided a test case in the comments section.

[Sean de Wolski]

How far off is it relatively?

I personally use this as a tester for MATLAB Coder generated MEX files:

% Unit test object
Tester = matlab.unittest.TestCase.forInteractiveUse;

% Assert within tolerance
assertEqual(Tester,foo(inputs),foo_mex(inputs),'RelTol',10^-15)

[Garren]

The relative differences are quite small, almost as small as 10^-10. For instance, for a known input that causes this issue, the norm(foo_mex_output) is usually around there.

However the relative percent change (which is what I am interested in) is equal to

   100*(foo_mex_output-foo_output)./(foo_output)

can result in values of over 100, which is far too great of a change in value for the current application of this MEX function.

[Garren]

As a test case, I've attached a .mat file containing five matrices. The contents are as follows and in this order:

in : The input matrix to the expm() (and the MEX version) for this specific case

U : The variable U, just before it is called by the division

V : The variable V, just before it is called by the division

F : The output after the division of the regular expm() function

mex_F : The output after the division of the MEX-version expm() function

Once the test file is downloaded, I've included the code for your convenience to extract the data.

   load('testCase.mat');
   [in, U, V, F, mex_F] = testCase{:};

The variables U and V are equivalent throughout both versions of this function. However, after this line is executed:

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

The values of the outputs (F and mex_F) are not equivalent and can be shown by '==' comparison. Not only that, but the percent change, which can be tested by running,

    100*(mex_F-F)./F

shows some elements have changed drastically. Some elements are inf which means that output values that should have been '0', were not.

Hopefully this will give a little more insight to my problem. Perhaps it has something to do with the properties of the input matrix, or could it just be a calculation error with the mldivide()?

I am running on Windows 7, and the compiler being used is Microsoft Windows SDK 7.1 (C++)

[John]

you could always use a Gaussian elimination function since I am pretty sure mldivide uses some sort of Gaussian elimination method in it.

[Matt J]

since I am pretty sure mldivide uses some sort of Gaussian elimination method in it.

No, it does not. See the flowcharts under Algorithm for Full Inputs and Algorithm for Sparse Inputs at

http://www.mathworks.com/help/matlab/ref/mldivide.html

[John]

Ahh thanks, that is helpful!

I had to code a replacement for mldivide for some code and I chose to use Gaussian Elimination and then back substitution. The Flow chart Matt provided should give you more directions towards what linear equation solvers would best fit your needs.

Answer 1

[Sean de Wolski]

The results are expected. If we look at the difference between F and Fmex it’s 2*eps(F)

eps(max(F(:)))
 ans =
 4.4409e-16

norm(F-mex_F)
 ans =
 9.9491e-16

Since some values in F are exactly 0:

any(F(:)==0)
 ans =
 1

Dividing by these will yield either nan (0/0) or +-inf (anything else/0).

The fair "percentage test" would be:

100*norm(F-mex_F)/norm(F)
 ans =
 3.1797e-14

Which is very small.

[Garren]

I understand, this makes much more sense. Thank you very much!

Answer 2

[Matt J]

It is generally not realistic to expect low percent error in two floating point calculations meant to result in zero or something close to zero. There are so many ways in floating point math that things meant to result in zero can cancel imperfectly, e.g.,

    >> x=[1.1.^(0:4) , -1.1.^(4:-1:0)]

    x =

        1.0000    1.1000    1.2100    1.3310    1.4641   -1.4641   -1.3310   -1.2100   -1.1000   -1.0000

    >> sum(x) %supposed to be zero

    ans =

       1.3323e-15

and obviously any such error relative to zero is going to be Inf percent.

[Garren]

I understand this, but my concern also lies with the fact that the results are exactly the same, up until it reaches the division point, and the previous lines of code involve a lot of matrix multiplication.

[Matt J]

All that means is that you did better than expected up until the "division point", not that you should expect as good later on.

Maybe you had integer matrix data which can be done in exact arithmetic if all you were doing is multiplying.