Inconsistency between Matlab eig() function and Matlab generated C code eig() function

Hello,

I'm looking to generate a C code from a Matlab function using Matlab in-built eig() function.

I have the following matrix :

M = [-0.0062    -0.1497    0.1435; ...
      0.8728    -0.6724    0.2994; ...
      1.7736    -0.4364   -0.0062];

The associated eigenvectors using Matlab in-built eig() function are :

V = [0.2184    0.3540   -0.0643; ...
     0.5088    0.4597   -0.7629; ...
     0.8327   -0.8145   -0.6433];

When using generated C code (currently using Simulink Accelerator) I get the following results :

Vc = [0.2019 - 0.0831i  -0.2295 - 0.2696i  -0.0572 + 0.0294i; ...
      0.4705 - 0.1937i  -0.2979 - 0.3500i  -0.6784 + 0.3489i; ...
      0.7701 - 0.3170i   0.5280 + 0.6202i  -0.5721 + 0.2943i];

Granted, the complex magnitude of the previous expression returns the same values as the absolute values of V, but I want to get the eigenvector associated to the minimal non-negative value of the expression 4*Vc(1,:).*Vc(3,:)-Vc(2,:).^2 and in that case, with Matlab, I end up with :

sol = [0.4685   -1.3646   -0.4164];

And with the generated C code, with :

solc = [0.3327 - 0.3298i   0.2179 - 1.3471i  -0.2422 + 0.3387i];

As you can see neither selecting the real or the imaginary part of the previous expression would allow for some consistency between Matlab code and generated C code. Taking the complexe magnitude of the previous expression would mean losing information about negative values.

Is there any workaround that I could use for both solutions to give consistent and comparable results ? Or do I need to scrap the idea of generating code while using the eig() function ?

Sorry in advance if there's some obvious mathematical solution to my issue, I'm not really comfortable working with eigenvalues/eigenvectors.

Answer 1

[Christine Tobler]

For code generation, the eig function doesn't have as many special-case treatments as the eig function in regular MATLAB. This means that both outputs are correct, the first one is u more convenient.

To see that both outputs are correct, you can compute M*V - V*D with D the diagonal matrix of eigenvalues. This is close to round-off error for both cases.

Now what are your options for dealing with this?

1) Generalize your output treatment to work for complex eigenvectors

While for this matrix all eigenvalues and eigenvectors are real, in general for a non-symmetric real matrix, it is possible that some of the eigenvalues and eigenvectors are complex.

Example:

rng(10);
[V, D] = eig(randn(3))

U = 
  -0.2079 - 0.6579i  -0.2079 + 0.6579i   0.4675 + 0.0000i
   0.6935 + 0.0000i   0.6935 + 0.0000i   0.3267 + 0.0000i
  -0.2042 + 0.0363i  -0.2042 - 0.0363i   0.8214 + 0.0000i

D = 
   0.2615 + 1.4204i   0.0000 + 0.0000i   0.0000 + 0.0000i
   0.0000 + 0.0000i   0.2615 - 1.4204i   0.0000 + 0.0000i
   0.0000 + 0.0000i   0.0000 + 0.0000i  -0.6271 + 0.0000i

So if you can find a way to make your formula on the rows of V also work for the case of complex eigenvectors, this will also make your algorithm safer when it's run in MATLAB without codegen.

What adjustment works will depend on where this formula is coming from, but at a guess it seems like 4*conj(Vc(1,:)).*Vc(3,:)-abs(Vc(2,:)).^2 matches up quite well with what is returned for V:

>> 4*conj(Vc(1,:)).*Vc(3,:)-abs(Vc(2,:)).^2
ans =
   0.4684 - 0.0000i  -1.3648 + 0.0001i  -0.4165 - 0.0001i
>> 4*conj(V(1,:)).*V(3,:)-abs(V(2,:)).^2
ans =
    0.4686   -1.3647   -0.4166

So my suggestion would be to use this formula instead and not worry about codegen having complex eigenvectors (after having verified that it makes sense in the context of what the formula represents, of course).

2) Include a LAPACK library reference for the generated code to call

EIG calls one of several methods of the LAPACK library depending on the input type. For code generation, we can't expect the common LAPACK implementations to compile on the targetted CPU, so a subset of these methods is implemented directly as part of MATLAB codegen (this is what was called in your case).

If you're on a CPU where a LAPACK library is available, you can add a coder.LAPACKCallback class to the code generation build so that the respective LAPACK function is called instead.

I haven't tried this myself, the linked doc pages suggest that possibly the LAPACK functions are only called for larger matrix sizes. So I would favor my first suggestion, but this is also an option.

[Hugoz]

Thank you for your suggestions. I will stick with the first one for now until I need to validate my algorithm.

For more background, I'm implementing an ellipse fitting method from : Halir, R. and Flusser, J. (1998). Numerically Stable Direct Least Squares Fitting Of Ellipses. with the algorithm presented on page 4.

For more consistency I'd rather make sure both the Matlab and the generated C code use the same libraries but it can wait.

Thanks again.

Edit :

@Christine Tobler : Do you mind clarifying where your guess of using 4*conj(Vc(1,:)).*Vc(3,:)-abs(Vc(2,:)).^2 is coming from ? Thanks.

[Christine Tobler]

My guess was that what's happening here is more-or-less an inner product where we have a vector of Vc being

[x; y; z]

and the result computed here is

[2*x; y].' * [2*z; y]

the previous formula would have worked assuming a, b, and c are real, but for complex we would want to use ctranspose (') instead of non-conjugate transpose (.'). My changes to the formula were making that adjustment.

That seems to match up with the comment in the paper that what this is computing is a'*C*a -> if a is complex, we would want to use a' (complex conjugate transpose) and not a.' (regular transpose).

[Hugoz]

Ok, that makes more sense now, thank you very much !

Answer 2

[Steven Lord]

From the documentation page for the eig function, specifically the C/C++ Code Generation item in the Extended Capabilities section, two items:

• V might represent a different basis of eigenvectors. This representation means that the eigenvector calculated by the generated code might be different in C and C++ code than in MATLAB. The eigenvalues in D might not be in the same order as in MATLAB. You can verify the V and D values by using the eigenvalue problem equation A*V = V*D.
• Outputs are complex.

There is no guarantee that the output from the C++ code will be identical to the output from MATLAB. Remember that if V is an eigenvector of A with eigenvalue d, then so is k*V for a non-zero scalar k.

M = [-0.0062    -0.1497    0.1435; ...
      0.8728    -0.6724    0.2994; ...
      1.7736    -0.4364   -0.0062];
[V, D] = eig(M)

V = 3×3

0.2184 0.3540 -0.0643 0.5088 0.4597 -0.7629 0.8327 -0.8145 -0.6433

D = 3×3

0.1923 0 0 0 -0.5308 0 0 0 -0.3463

Vc = [0.2019 - 0.0831i  -0.2295 - 0.2696i  -0.0572 + 0.0294i; ...
      0.4705 - 0.1937i  -0.2979 - 0.3500i  -0.6784 + 0.3489i; ...
      0.7701 - 0.3170i   0.5280 + 0.6202i  -0.5721 + 0.2943i];
checkV = M*V-V*D

checkV = 3×3

1.0e-15 * 0.0208 0 -0.1145 0.0694 0.1110 -0.2776 0.0555 -0.7216 -0.4996

checkVc = M*Vc-Vc*D

checkVc = 
   1.0e-04 *

   0.0717 - 0.0089i  -0.2251 - 0.2723i   0.0652 + 0.0060i
  -0.3454 + 0.4445i  -0.2935 - 0.4442i   0.1670 - 0.0325i
  -0.6757 + 0.5517i  -0.7084 - 0.9148i   0.3378 - 0.2528i

It may seem like checkVc is not that close to 0, but if you'd given us more than 4 decimal places of the output from the C++ code I'd bet it would contain values that are similarly close to 0 as the ones in checkV. If we use the 4 decimal place approximation from your post instead of the full double precision V from the computation:

V2 = [0.2184    0.3540   -0.0643; ...
      0.5088    0.4597   -0.7629; ...
      0.8327   -0.8145   -0.6433];
howCloseAreVandV2 = V-V2

howCloseAreVandV2 = 3×3

1.0e-04 * -0.3498 -0.0532 -0.4377 -0.2747 -0.4187 0.2638 0.4680 0.0496 -0.3795

checkV2 = M*V2-V2*D

checkV2 = 3×3

1.0e-03 * -0.0178 -0.0042 0.0243 -0.0072 -0.0028 0.0582 0.0593 -0.0114 0.1020

Remember how above I said that eigenvectors aren't unique? Let's divide each element in Vc by its corresponding element in V and see what values we get for k.

k = Vc./V

k = 
   0.9246 - 0.3806i  -0.6483 - 0.7616i   0.8890 - 0.4569i
   0.9248 - 0.3807i  -0.6481 - 0.7614i   0.8893 - 0.4573i
   0.9248 - 0.3807i  -0.6483 - 0.7615i   0.8893 - 0.4575i

So if we scale the first eigenvector does it still satisfy the eigenvalue/eigenvector relationship?

KV = k(1, 1)*V(:, 1)

KV = 
   0.2019 - 0.0831i
   0.4704 - 0.1936i
   0.7700 - 0.3169i

checkKV = M*KV-KV*D(1, 1)

checkKV = 
   1.0e-16 *

   0.2082 - 0.0347i
   0.6939 - 0.3469i
   0.2776 - 0.1388i

Looks like it.

[Hugoz]

Thanks for you detailed answer. Since all related questions highlights differences with either the sign of the eigenvectors, their order, or the presence of a zero complex part, I was indeed surprised to get a non-zero complex part from the generated C code.

I don't doubt however that those two basis are both correct. My issue is how to go about comparing the two sets of eigenvectors to my condition in a consistent manner.

Answer 3

[William Rose]

@Hugoz, I am not sure how you generated the C code, but I suspect the C code is wrong. Matlab's eig() is extremely well tested and verified by real world use. Maybe you can go to a C site for assistance debugging the C code.

[Steven Lord]

Can it give different answers? Yes.

Can it give incorrect answers? If it does that would be a bug, but I don't see any evidence that its answers are incorrect.

[Hugoz]

I used Simulink Accelerator (which generate a C code) and extracted the results to the workspace. I have yet to generate the standalone C code with Matlab Coder + Embedded Coder. As soon as it's done, compiled and tested, I'll let you know if the results are different.