Both eig() and eigs() function calculate different eigenvalues depending on optional input values

The goal of my problem is to calculate the eigenfrequencies of a dynamic system by solving the generalized eigenvalue problem with .

Afterwards the eigenfrequencies in Hz are obtained by calculating .

The reference eigenfrequencies were calculated by two different finite element programs.

My observations are as follows:

• The eig() function calculates wrong eigenvalues with the default algorithm (Cholesky) and the right eigenvalues with the QZ-algorithm (Due to the condition?)
• The eigs() function calculates the right eigenvalues if the number of requested eigenvalues is and the wrong eigenvalues if
• Both function calculate the same wrong eigenvalues
• Only the first six eigenvalues are wrong, the other eigenvalues are very close to the right values

How can you explain this behaviour? If it is due to the condition of M, is there a threshold for the condition number until the default eig() will work?

Why does the number of requested eigenvalues in eigs() have such an impact on the result?

A MWE is shown below:

% script shows an MWE that the eig() and eigs() functions aren't very
% robust to badly conditioned matrices
% the reference eigenfrequencies are validated by two different FE
% programs:
% [9.097; 34.864; 100.236; 138.499; 206.162; 261.210; 689.248; 2357.613; 
% 4109.622; 9311.9419; 9455.152; 10254.968]

clear

%% load system matrices

load('matricesMK.mat') 

%% Solve an eigenvalue analysis
% the eig() function calculates the right results with the 'qz' algorithm.

[Veig,Deig] = eig(full(K),full(M),'qz');
[ome_eig,tmp] = sort(sqrt(diag(Deig)),'ascend');
frq_eig_qz = ome_eig/2/pi;

% By default the 'chol' algorithm is used and the results are wrong

[Veig,Deig] = eig(full(K),full(M));
[ome_eig,tmp] = sort(sqrt(diag(Deig)),'ascend');
frq_eig_chol = ome_eig/2/pi;

% the eigs() function calculates right eigenvalues for a number of
% eigenvalues k<18
k=17;
[V,D] = eigs(K,M,k,'sm');
[ome,tmp] = sort(sqrt(diag(D)),'ascend');
frq_eigs_right = ome/2/pi;

% the eigs() function calculates right eigenvalues for a number of
% eigenvalues k>=18
k=18;
[V,D] = eigs(K,M,k,'sm');
[ome,tmp] = sort(sqrt(diag(D)),'ascend');
frq_eigs_wrong = ome/2/pi;

% Further observations:
% Only the first six eigenvalues are wrong, the other eigenvalues are very
% close to the right values

[Bruno Luong]

You are dealing with numerically singular matrices, so you shouldn't' expect any consistent results with finite precision calculation.

>> diag(Deig)

ans =

   1.0e+21 *

  -0.000000000000000
   0.000000000000000
   0.000000000000000
   0.000000000000001
   0.000000000000001
   0.000000000000002
   0.000000000000019
   0.000000000000219
   0.000000000000667
   0.000000000003423
   0.000000000003529
   0.000000000004152
   0.000000000008791
   0.000000000009242
   0.000000000011907
   0.000000000018025
   0.000000001491550
   0.000000002833869
   0.000000003393521
   0.000000003890927
   0.000000004247819
   0.000000006339370
   0.000000007487080
   0.000000007880103
   0.000000007988986
   0.000000008324769
   0.000000008330993
   0.000000008331892
   0.000000008334342
   0.000000011109544
   0.416504353653501
   0.416504353653522
   0.422208709630153
   0.440299455676799
   0.498737193433392
   1.223316108193663


>> cond(K*inv(M))
Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND =  1.643016e-17. 
 
Warning: Using CONDEST instead of COND for sparse matrix. 
> In cond (line 25) 

ans =

     1.725316053089528e+19

[Johannes R.]

Thank you for your answer @Bruno Luong!

I know that the matrix M is badly conditioned but I still can't fully understand the behaviour.

The condition number is pretty high but still not Inf. Do we need to consider the matrix as singular?

The algorithms seem to converge in both cases and they converge to the right values by simply making small changes.

Why does the number of requested eigenvalues in eigs() change to results? And why are the results the same wrong results as for eig()? Is the algorithm switching to 'chol' at this point?

[Bruno Luong]

When your matrix is ill-conditioned, any linear algebra algorithm is likely to be affected by round-off error.

Any matrix with condition number > 1e6-1e7 would be considered as singular numerically.

Please google about "condition numbers", there are a lot of writting about the topic.

[Johannes R.]

Alright, I didn't expect it to be so sensitive with respect to the condition.

Thank you for your help!

Answer 1

[Christine Tobler]

The problems seen here are due to M being symmetric positive semi-definite. Such a matrix is numerically singular, which can have quite a bit of unfortunate effects.

EIG - 'chol' vs 'qz' and positive semi-definite M

In general, the 'qz' algorithm is more robust than the 'chol' algorithm, it works even if K and/or M are singular.

However, for symmetric problems (where K is symmetric and M is symmetric positive definite), EIG makes the promise that the outputs will have real eigenvalues and M-orthogonal eigenvectors (V'*M*V = I). This is not guaranteed for the 'qz' algorithm, which works for any non-symmetric system.

This is where the 'chol' option comes in: The Cholesky decomposition of M is computed, and this is used to rewrite the generalized eigenvalue problem as a simple, symmetric eigenvalue problem:

R = chol(full(M));  % M = R'*R
RKR = R'\K/R;       % K * x = lambda * (R'*R) * x   <=>   R'\ K / R * y = lambda * y, with y = R*x.
[Veig,Deig] = eig((RKR + RKR')/2); % Make sure input is exactly symmetric
% Note: Veig would need to be adjusted, too.
[ome_eig,tmp] = sort(sqrt(diag(Deig)),'ascend');
frq_eig_chol_by_hand = ome_eig/2/pi;

You'll see that this is quite similar to what EIG computes with the 'chol' option.

In principle, CHOL will error unless M is symmetric positive definite. However, for a symmetric positive semi-definite M, chances are 50-50 on whether it sees this as a badly conditioned sps-d matrix, or rejects it as an indefinite matrix (in which case it would fall back to 'qz').

In short, the problem is that the 'chol' option solves a linear system with R = chol(M), which introduces large errors when M is ill-conditioned.

EIGS - M-orthogonal Krylov subspaces and cut-off choice

Now in EIGS(K, M, k, 'sm'), the algorithm used doesn't require inverting M: It only involves backslash with K and matrix multiply with M. What's done internally is that an M-orthogonal basis of a Krylov subspace is constructed, then K is projected into that basis and EIG is applied to this smaller matrix. This is safe to do for symmetric positive semi-definite matrices.

So why does it fail for k == 18, in such a similar way to how EIG failed? That's because it calls EIG(full(K), full(M)) in this case.

The subspace I just mentioned above has a size that's heuristically set to 2*k, which is usually significantly smaller than the size of the whole matrix (this is the main intention of EIGS, quickly computing few eigenvalues of a large matrix). When 2*k >= the matrix size, computing EIG is quicker and in almost all cases more accurate than doing the explicit algorithm.

You can avoid this issue at k == 18 by setting the 'SubspaceDimension' Name-Value pair to any value less than 36. Use the 'Display' Name-Value pair to see if EIGS is calling into EIG or starting a Krylov-Schur iteration.

[Walter Roberson]

You really know the linear algebra routines well, Christine !! When I need to understand the routines better, I refer to your posts!

[Johannes R.]

Thanks a lot for the detailed view into the matlab routines!

It really helped me to understand the shown behaviour and now I know where to focus on when handling such kind of problems.