Computational Performance and Accuracy Comparison of Standard Eigenvalue Problem eig(A) and Generalized Eigenvalue Problem eig(A,B)

I have been thinking about generalized eig(A) and standard eigenvalue eig (A,B) problem. I did a finite element convergence study in MATLAB using both these two methods by increasing the number of elements upto 1500. I compared both method in MATLAB for a shaft system with gyroscopic effect and without gyroscopic effect (0 speed). I have found that standard eigenvalue problem is giving exactly same result in state space form for 0 zero speed form compared with eig(K,M). However, generalized eigenvalue problem can give randomly wrong modes after 360 number of element. On the other hand, standard eigenvalue problem is so faster than generalized eigenvalue problem. For 1200 elements, generalized takes 10 hours whereas standard takes 42 minutes.

Does anyone have an idea about this phenomenon? How can MATLAB give different results when using standard and generalized eigenvalue? I found that using standard eigenvalue problem gives better result than generalized eigenvalue problem in terms of computational efficiency and accuracy.

[Christine Tobler]

I did not understand if your problem is a standard or a generalized eigenvalue problem? Giving different results should be expected, as these are two different problems. It's also expected for generalized to be slower than standard, as this is a harder problem.

Could you describe how you are determining that one of the results is wrong?

Also, the times you are describing seem very slow: On my machine, a standard eigenvalue problem of size 1200-by-1200 takes 2 seconds, while the related generalized problem takes 23 seconds.

 >> A = randn(1200);
 >> tic; [U, D] = eig(A); toc
 Elapsed time is 1.927305 seconds.

 >> A = randn(1200); B = randn(1200);
 >> tic; [U, D] = eig(A, B); toc
 Elapsed time is 23.093472 seconds.

If possible, could you upload the matrices you are using?

[Ali Tatar]

I can formulate quadratic eigenvalue problem using generalized and standard form in state space representation as you see below;

A1=[(C+G) M; M zeros(sdof)]; % Generalized Eigenvalue Problem

B1=[K zeros(sdof); zeros(sdof) -M];

[eigvec_right,eigval,eigvec_left]=eig(B1,-A1);

AA=[zeros(sdof) eye(sdof); -M\K -M\(G+C)]; % Standard Eigenvalue Problem

[eigvec_right,eigval,eigvec_left]=eig(AA);

where C is symmetric damping matrix, G is skew symmetric gyroscopic matrix, M is symmetric mass matrix and K is symmetric stiffness matrix.

The size of the matrix is becoming 14412 if 1200 beam elements is used. The formula is size=(nel+1)*ndof*2

where nel is number of elements, ndof is number of degrees of freedom. nel is 1200 and ndof is 6. Due to using state space representation size is becoming 2n*2n.

Theoretically, these two methods should give the same result. However, generalized eigenvalue problem computes randomly new modes at higher number of elements, in general after 360 elements. For this reason, I have been using standard eigenvalue problem rather generalized one. But I wonder why generalized eigenvalue problem gives wrong result and is slower than standard eigenvalue problem.

[Christine Tobler]

For the generalized problem, you are computing Ug, Dg to solve

 B1 * Ug = - A1 * Ug * Dg

For the standard problem, you are computing Us, Ds to solve

 AA * Us = Us * Ds

From these two equations, I would assume that AA = -A1 \ B1, but this does not match your definition of AA. It looks more like you are setting AA = B1 \ A1, but that also doesn't match the definition, since this would result in

 AA = [K \ (C + G), K \ M; -eye(sdof), zeros(sdof)];

It would help if you could upload an example script that constructs the matrix and shows the problem. Some more questions about the problem: Are M, K symmetric positive (or negative) definite? Are any of the matrices singular?

[Ali Tatar]

These two type of matrix construction can be found in Daniel Inmann Engineering Vibration Book and Michael Friswell`s Dynamics of Rotating Machine Book. I can say that they are identical to each other. Also, I am getting the same result with two methods for less than 360 number of elements (14412 matrix size). Therefore, there is no hesitation for the state space construction. The problem is here obvious computation time difference between these two methods and getting wrong result in generalized one for large eigenvalue problems. M, K and C matrices are symmetric. G matrix is skew symmetric. There is no singular matrix.

Is it possible to learn which eigensolver and numerical algorithm is used in generalized and standard eigenvalue problem in MATLAB?

[Christine Tobler]

I'm not sure if we are talking about the same thing, in terms of the answers between the two systems having to match. The function EIG will return a solution of the (generalized) eigenvalue problem; you can check if the solution is correct by computing the residual.

Here's how to compute the residuals for both cases. Could you try out the residuals for the case where standard and generalized were not agreeing?

 % Standard eigenproblem:
 [U, D] = eig(A);
 residual = A*U - U*D;
 absRes = norm(residual)
 relRes = norm(residual / D)

 % Generalized eigenproblem:
 [U, D] = eig(A, B);
 residual = A*U - B*U*D;
 absRes = norm(residual)
 relRes = norm(residual / D)

About the algorithms used in EIG: MATLAB uses LAPACK functions for linear algebra:

- For the standard problem, there is a specialized algorithm for the symmetric case (only used if issymmetric(A) is true, not if the matrix is slightly off). Some of the algorithms are described here: http://www.netlib.org/lapack/lug/node48.html, I'm not sure off-hand which of them is used in EIG.

- For the generalized problem (eig(A, B)) a specialized algorithm is used if B is symmetric and positive definite. If B is symmetric but not positive definite, the non-symmetric variant is used: the QZ algorithm. From the run-time of your call, I would assume that your second input to eig was not symmetric positive definite. Maybe this could be fixed by adjusting the sign, or switching the order of inputs (using eig(A1, B1) instead of eig(B1, A1), and then transforming the eigenvalues as needed).

Mathematically, eig(A, B) is equivalent to eig(B\A). However, the computation of B\A introduces numerical error, so eig(A, B) is typically more numerically accurate.

[Ali Tatar]

Generalized Eigenvalue Solution:

Nel 30

absRes =

    0.0024

relRes =

   6.8796e-06

Nel 60

absRes =

    0.0087

relRes =

   3.6068e-05

Nel 360

relRes =

    0.0014

absRes =

    0.3648

Standard Eigenvalue Solution

Nel 30

absRes =

   1.4700e-06

relRes =

   3.0110e-08

Nel 60

absRes =

   4.9571e-06

relRes =

   1.0154e-07

Nel 360

absRes =

   1.3813e-04

relRes =

   2.8286e-06

The problem is here that when I increase the number of beam elements in my finite element code, generalized eigenvalue solver gives wrong result and is creating new eigenvalues at 360 beam element.(nel = number of element). It is not converging.

[Christine Tobler]

From what I can see here, this could be a numerical issue, or it could be an effect of a generalized eigenvalue problem that is badly conditioned. It is hard to say more without seeing the matrices and outputs.

If you are able to send us a MATLAB script that generates these matrices, or a .mat file containing the matrices, we could take a closer look. In this case, please contact MathWorks tech support.

[Ali Tatar]

By the way, I checked the absRes, relRes values of generalized eigenvalue solution for 300 number of elements, which gives the right result and closest value to the 360 number of elements. Something is happening when it becomes large eigenvalue problem.

Nel 300 Standard Eigenvalue Solution

absRes =

    0.2527

relRes =

   4.8100e-04

In my opinion it can be a numerical issue of solution routine. Unfortunately, I am not able to send the matlab script. I can say that it is a finite element code for a beam. It automatically divides whole beam into the elements. If you have a finite element code, you can easily run this analysis using generalized eigenvalue solution.

[Christine Tobler]

Sorry I couldn't help you more, there is not much to be done when operating blindly. While an issue with the EIG function is not impossible, I think it's much more likely that the input to EIG is badly conditioned.

A possible problem with finite element code can be if boundary conditions are not removed in the correct way: This will make the problem eig(A, B) have no solution - but with numeric errors, it can make it just practically unsolvable, and EIG will try its best to return something.