EIG incorrect for large matrices when compiled into C on speedgoat target computer
3 views (last 30 days)
Show older comments
Patrick Franks on 19 Dec 2018
Commented: Delaney 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?
Delaney on 20 Jul 2019
Edited: Delaney 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)
Delaney on 21 Jul 2019
Edited: Delaney on 21 Jul 2019
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:
[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);
[-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];
[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);
[ 0.0406, 0.0451, -0.9982;
0.6582, -0.7528, -0.0072;
-0.7518, -0.6567, -0.0602];
[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
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.
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.
Delaney on 22 Jul 2019
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):
[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 on 1 Aug 2019
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.
Find more on Linear Algebra 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!