MATLAB Answers

0

EIG incorrect for large matrices when compiled into C on speedgoat target computer

Asked by Patrick Franks on 19 Dec 2018
Latest activity Commented on by Delaney Miller on 1 Aug 2019
I have been using Matlab R2017b to create simulink files that are compiled onto a speedgoat target computer. In the simulink file, in a matlab function block, we use the function eig to find the eigenvectors and eigenvalues of a matrix A. When A is 4 by 4 or smaller, this works fine.
When A is larger than 4x4 (i.e. 5x5, 6x6), eig does not correctly return the eigenvalues of the matrix. That is, for the same matrix A, the eigenvalues generated by eig on the target computer are different than those generated on the target computer (i.e. calling eig in the workspace normally).
Is this related to the code being compiled into C/C++? Has anyone experienced this issue?
Thanks,
Patrick

  4 Comments

Show 1 older comment
Can you please provide some more details (ideally reproduction steps, and the outputs of EIG both in MATLAB and after being compiled from C++)?
The difference in eigenvalues is on the order of 10^26.
Often in C++ it will return NaN's or extremely large numbers. In MATLAB the outputs are correct.
The generation of the matrix is a little messy, but the matrix is a 5x5 positive definite symmetric matrix A. Here's an example of the code and output.
[B,D] = eig(A).
In Matlab, B(1) = -0.361.
In C++, B(1) = NaN.
I suspect that the C++ version of eig is not getting the correct A as input. Maybe dump A to a file (with the C++ implementation) and verify that it has the intended values.

Sign in to comment.

3 Answers

Answer by Delaney Miller on 20 Jul 2019
Edited by Walter Roberson
on 21 Jul 2019

[B, D] = eig(C, eye(N), 'qz'); % N = dimension of matrix
for i = 1:N
B(:,i) = B(:,i) ./ norm(B(:,i)); % normalize eigenvectors (columns)
end
;)

  2 Comments

Can you provide a small code snippet that we can run to create C and define N or edit your answer and attach a MAT-file containing C and N to it? That way we can try to reproduce the behavior you described.
A section of the code requires the eigendecomposition of a covariance matrix C that is real and symmetric. Here is the block of code (in my case, C is a 3x3 matrix):
N = 3;
C = triu(C) + triu(C,1)'; % enforce symmetry
[Btemp,Dtemp] = eig(C); % eigen decomposition (Btemp is normalized)
D = sqrt(diag(Dtemp)); % D is a vector of standard deviations now
invsqrtC = real(Btemp * diag(D.^-1) * Btemp'); % sqrt(C^-1)
D = real(D);
B = real(Btemp);
Previously, when C was a 4x4 matrix, the Btemp and Dtemp matrices computed by eig(C) had real entries (version running was converted C/C++ code). However, I noticed when I changed C to a 3x3 matrix (all dimensions scaled accordingly), the eig(C) function produced complex eigenvectors in the Btemp matrix, and the eigenvalues in the Dtemp matrix did not agree with the same code computed in Matlab directly. Here is an example:
C =
[0.9736, 0.0001, 0.0066;
0.0001, 0.8655, 0.0115;
0.0066, 0.0115, 0.8628];
On the Speedgoat target computer (C/C++ code version), the following Btemp and Dtemp matrices are computed:
% C/C++ code version on Speedgoat target computer
[Btemp,Dtemp] = eig(C);
Btemp =
[-0.5415 + 0.0000i, 0.0121 + 0.2506i, -0.6582 - 0.0046i;
-0.5463 + 0.0000i, 0.8889 - 0.1470i, 0.3861 - 0.3385i;
-0.8785 + 0.0000i, -0.0878 - 0.1667i, 0.4377 + 0.0334i];
Dtemp =
[0.6670, 0, 0;
0, 0.6590, 0;
0, 0, 0.6590];
When the same calculation is performed in Matlab with the same C matrix, the following Btemp and Dtemp matrices are computed:
% Matlab version on host computer
[Btemp,Dtemp] = eig(C);
Btemp =
[ 0.0406, 0.0451, -0.9982;
0.6582, -0.7528, -0.0072;
-0.7518, -0.6567, -0.0602];
Dtemp =
[0.8524, 0, 0;
0, 0.8755, 0;
0, 0, 0.9740];
I would like to note that with a 4x4 C matrix this was not a problem. The eigenvalues and eigenvectors appeared in a different order in the Speedgoat (C/C++ code) version and the Matlab version, but all of the same eigenvalues were computed. I solved this problem myself by replacing the following line:
[Btemp, Dtemp] = eig(C);
With a version that ensures the 'qz' eigendecomposition algorithm:
N = 3;
[Btemp, Dtemp = eig(C, eye(N), 'qz');
for i = 1:N
Btemp(:,i) = Btemp(:,i)./norm(Btemp(:,i)); % normalize eigenvectors
end
When I do this, I get the same results on the Speedgoat target computer (C code) as the Matlab version I calculate with eig(C) (except for a sign reversal on the eigenvectors, but that's not important). I believe there may be an issue with how the 'schur' algorithm is converted into C/C++ code using Matlab Coder, specifically with odd sized matrices. In the original question, both 5x5 and 9x9 matrices were tried. However, I've confirmed that 4x4 matrices do not experience ths problem.

Sign in to comment.


Answer by Mike Hosea
on 22 Jul 2019

Bit difficult to debug from here. Is your code calling the initialize and terminate functions for the library you are generating?
As for complex eigenvectors, for Hermitian A we have been generating code that performs [V,D] = schur(A,'complex') and then zeros the off-diagonal elements of D. That is to say, code generation has been treating every problem as a complex one in EIG, which in the general case it must because complexity is immutable in the generated C code and depends on input data values. That might explain the occasional appearance of complex eigenvectors. We're planning to change this to do [V,D] = schur(A,'real') instead when A is real and symmetric. Since you know your input matrix is real and symmetric, you could just do that yourself instead of calling EIG.
Of course none of that explains incorrect results. I'm not aware of anything that would cause norm(A*V - V*D) to be larger than expected, and I'm not seeing any memory overruns or anything like that when debugging generated code that doesn't call into a LAPACK library, such as would be the case for a stand-alone target.
We do regularly find bugs in C compilers around here. It is worth generating code with different optimization levels to see what happens.
Finally, I believe we supported calling LAPACK on a stand-alone target in 17b. If you can get a LAPACK library built for your target, that might be a workaround and and also deliver performance benefits on larger problems. I'm not personally familiar with the process, but I know it is documented.

  0 Comments

Sign in to comment.


Answer by Delaney Miller on 22 Jul 2019

Hi Mike,
Thanks for the help. What do you mean by your first question? We are not initializing/terminating the eig function as we are using it within a larger function (Simulink fcn block). We are using a Speedgoat target machine (https://www.speedgoat.com/learn-support/simulink-real-time-workflow) and I believe our driver blocks are working correctly as it isn't throwing any errors. Maybe this is an error on Speedgoat's end when the C code is compiled.
I'm aware of the default settings with 'schur' to the complex eigendecomposition. When I use schur(C,'real'), I find that I similarly get a completely wrong eigendecomposition. It's just all real instead of complex and the Dtemp matrix (of eigenvalues) has non-zero entries for non-diagonal indices. Here's an example (slightly different C matrix):
Dtemp =
[6.9407 -32.0045 -0.9886;
57.0459 6.9407 2.4762;
0 0 6.9407];
So yeah, it seems the issue is with the 'schur' algorithm as 'qz' works fine. Because we only run this computation every 15 minutes or so within the context of a larger biomechanics data collection, it doesn't necessarily make sense for us to go to the trouble of building a LAPACK library. I understand that 'qz' algorithm is not as efficient (esp. for a real symmetric / Hermitian matrix) but it seems to be accurate so we'll stick with that solution for the time being.
Appreciate you looking into this,
Delaney

  4 Comments

Show 1 older comment
BTW, I've created an internal bug report to investigate this as best we can. I doubt the problem is actually with the SCHUR algorithm itself, rather either with a transform being done when generating code from that algorithm or, as I suggested before, a C compiler bug. Hopefully we'll be able to sort it out.
I haven't been able to reproduce the issue on desktop. However, there has been some discussion here about whether the fastmath optimization ( /fp:fast ) might be playing a role. Apparently /fp:fast is used for SpeedGoat targets. You might try deleting that from the generated .mk file and rebuilding to see if it changes anything.
Makes sense that the issue could be with Speedgoat. Right now using 'qz' seems to be fine, but if I have the chance I'll try deleting the fastmath optimization and let you know what happens.

Sign in to comment.