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

17 次查看（过去 30 天）
Hugoz
Hugoz 2024-1-19
评论： Hugoz 2024-1-23
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.

请先登录，再回答此问题。

采纳的回答

Christine Tobler
Christine Tobler 2024-1-19
编辑：Christine Tobler 2024-1-19
Ran in:
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 个评论
显示 1更早的评论
Christine Tobler
Christine Tobler 2024-1-23
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).

请先登录，再进行评论。

更多回答（2 个）

Steven Lord
Steven Lord 2024-1-19
Ran in:
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 个评论
Hugoz
Hugoz 2024-1-19
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.

请先登录，再进行评论。

William Rose
William Rose 2024-1-19
@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 个评论
Steven Lord
Steven Lord 2024-1-19
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 2024-1-19
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.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Mathematics and Optimization 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by