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

27 views (last 30 days)
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.

Accepted Answer

Christine Tobler
Christine Tobler on 19 Jan 2024
Edited: Christine Tobler on 19 Jan 2024
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.
  3 Comments
Christine Tobler
Christine Tobler on 23 Jan 2024
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).

Sign in to comment.

More Answers (2)

Steven Lord
Steven Lord on 19 Jan 2024
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.
  1 Comment
Hugoz
Hugoz on 19 Jan 2024
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.

Sign in to comment.


William Rose
William Rose on 19 Jan 2024
@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.
  2 Comments
Steven Lord
Steven Lord on 19 Jan 2024
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
Hugoz on 19 Jan 2024
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.

Sign in to comment.

Categories

Find more on Mathematics and Optimization in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!