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 
.
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
4 个评论
Bruno Luong
2021-3-17
编辑：Bruno Luong
2021-3-17
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.
2021-3-17
Bruno Luong
2021-3-17
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.
2021-3-18
采纳的回答
Christine Tobler
2021-3-24
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.
2 个评论
Walter Roberson
2021-3-24
编辑：Walter Roberson
2021-3-24
You really know the linear algebra routines well, Christine !! When I need to understand the routines better, I refer to your posts!
Johannes R.
2021-3-25
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Linear Algebra 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
