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

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

[Walter Roberson]

how much difference ? is this rounding error or is it very different ?

[Christine Tobler]

Can you please provide some more details (ideally reproduction steps, and the outputs of EIG both in MATLAB and after being compiled from C++)?

[Patrick Franks]

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.

[Matt J]

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.

Answer 1

[Delaney]

[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

;)

[Steven Lord]

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.

[Delaney]

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.

Answer 2

[Mike Hosea]

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.

Answer 3

[Delaney]

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

[Mike Hosea]

I missed that you were working in Simulink. When generating a library for deployment on some other target with MATLAB Coder, you're required to call an initialize function before you actually call one of the (other) entry point functions of the library. Failing to do that can cause some very confusing situations. That's what I was referring to, but it's not relevant here.

For a long time all we had was QZ. IIRC, the principal motivation for using SCHUR was to ensure that the eigenvectors would be orthogonal when the input is Hermitian. The code snippet above assumes it, for example. I don't remember the details, but we had cases where they weren't, IIRC because of repeated eigenvalues. With QZ you should probably use /Btemp instead of *Btemp' to be on the safe side.

[Mike Hosea]

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.

[Mike Hosea]

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.

[Delaney]

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.