When Should Generalized eig() Solver be Used for Ordinary eig() Problems?

61 visualizzazioni (ultimi 30 giorni)
Paul
Paul il 13 Apr 2024 alle 21:35
Commentato: Bruno Luong il 14 Apr 2024 alle 20:12
Ran in:
Consider the following matrix:
z = sqrt(sym(2));
z1 = sqrt(sym(3));
Asym = [-z, -z^2/z1; z1, z]
Asym = 
It has a repeated eigenvalue at the origin
eig(Asym)
ans = 
Enter the matrix as a double
format short e
z = sqrt((2));
z1 = sqrt((3));
Anum = [-z, -z^2/z1; z1, z]
Anum = 2x2
-1.4142e+00 -1.1547e+00 1.7321e+00 1.4142e+00
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
It is nilpotent in double precision
Anum*Anum
ans = 2x2
0 0 0 0
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
Its eigenvalues using ordinary eig() are on the order of 1e-8.
eig(Anum)
ans =
2.5025e-17 + 2.0134e-08i 2.5025e-17 - 2.0134e-08i
Its eigenvalues solving as a generalized eigenvalue problem (suggested by Torsten)
eig(Anum,eye(2))
ans = 2x1
1.0e+00 * -3.3920e-17 -7.0217e-17
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
Why is using the generalized eig so much more accurate in this case? More generally, are there guidelines on when to use the generalized eig solver for ordinary eigenvalue problems based on characteristics of A?
  3 Commenti
Mostra 1 commento meno recente
Torsten
Torsten il 13 Apr 2024 alle 23:32
Modificato: Torsten il 13 Apr 2024 alle 23:33
This matrix do not have two eigenvalues (or one with multiplicity 2). It has just one (0) with multiplicity 1, and the matrix is not diagonalisable.
Better geometric multiplicity - the algebraic multiplicity is 2 since the characteristic polynomial is lambda^2.
Bruno Luong
Bruno Luong il 14 Apr 2024 alle 5:46
Modificato: Bruno Luong il 14 Apr 2024 alle 8:41
Ran in:
It you call EIG with eigen vector out they are the same (with some numerical error). They do not form a basis of R^2. So eigen values 1 and 2 ireturned by EIG is the "same" eigen value.correspond to the same 1D sibspace span by this vector,
format short e
z = sqrt((2));
z1 = sqrt((3));
Anum = [-z, -z^2/z1; z1, z];
[V,Dcall1]=eig(Anum,'vector')
V =
-6.3246e-01 + 9.0044e-09i -6.3246e-01 - 9.0044e-09i 7.7460e-01 + 0.0000e+00i 7.7460e-01 + 0.0000e+00i
Dcall1 =
2.5025e-17 + 2.0134e-08i 2.5025e-17 - 2.0134e-08i
Dcall2=eig(Anum)
Dcall2 =
2.5025e-17 + 2.0134e-08i 2.5025e-17 - 2.0134e-08i
rank(V,1E-7)
ans =
1
isequal(Dcall1,Dcall2)
ans = logical
1
Only the symoclic EIG returns correctly a single eigen vector:
z = sqrt(sym(2));
z1 = sqrt(sym(3));
Asym = [-z, -z^2/z1; z1, z]
Asym = 
[V,D,p]=eig(Asym)
V = 
D = 
p =
1

Accedi per commentare.

Accedi per rispondere a questa domanda.

Risposta accettata

Bruno Luong
Bruno Luong il 14 Apr 2024 alle 8:17
Modificato: Bruno Luong il 14 Apr 2024 alle 8:51
Ran in:
Actually I think eig with 2 inputs is more approprite for multiple order eigen value and the degenerate case to non trivial jordan form The Schur decomposition (using diferent left and right orthogonal on left and right side of the decompositopns) reveals the jordan form:
z = sqrt((2));
z1 = sqrt((3));
A = [-z, -z^2/z1; z1, z];
I = eye(2);
[TA,TI,Q,Z] = qz(A,I);
norm(TI - I) % TI is pratically unchange from I
ans = 2.2208e-16
TA
TA =
0.0000 - 0.0000i -2.8868 - 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i
d = diag(TA);
errd = norm(d,Inf) % close to 0, the eigen value of A
errd = 3.0417e-08
TA1termedia2 = TA(1,2) % This is the off diagonal term of the Jordan form
TA1termedia2 = -2.8868 - 0.0000i
Now how on earth by which miracle an intermediate errd of 3e-8 lead to even more accurate (the error get squared) final eigen values
eig(A,I,"qz")
ans = 2x1
1.0e-16 * -0.3392 -0.7022
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
I have no idea !
QR algorithm used by EIG with single input can be viewed as some sort of power iterative method and less performant to reveal jordan block structure of the matrix.
Not tested but note that another algorithm branch of generalized eig is possible using cholesky decompolizion in case both matrix is symmetric and the second one is positive, I think it is also robust to multiplicity as the theory says the eigen vectors are orthogonal and numerical method will enforce it.
So I change my opinion : better use EIG with 2 arguments (generalized) for repeat eigen value case. Torsen have found a good trick here.
  3 Commenti
Mostra 1 commento meno recente
Bruno Luong
Bruno Luong il 14 Apr 2024 alle 16:18
Modificato: Bruno Luong il 14 Apr 2024 alle 18:36
Do you think it more appropriate only if the geometric multiplicity is less than the algebraic multiplicity, as in this particular problem? Or more appropriate even in cases where the geometric multplicity equals algebraic multiplicity?
The example in this thread seems to indicate it is more robust even for eigenvalue with 'geometrical) multiplicity > 1.
Would it matter if A has other, non-repeated eigenvalues?
Normally I would say no, distinct eigen values has "local" effect only. I could have sol ecross talk if the resêctive eigen spaces has very small angles.
Is there a way to know when to use the trick? That is, is there some characteristic(s) of A that can be robustly assessed to determine, a priori, that eig(A,I) is preferred over eig(A)?
I perdonally don't know any way.
Alternatively, is there any reason to ever not use eig(A,I), assuming that accuracy of the result is the only consideration?
Computation time? (oh but you say only accuracy matter). Actually I'm not sure if there is a conscensus about accuracy in general.
Can eig(A,I) ever give a less accurate result than eig(A)?
Someone should show an example then I could give some sriou s though. Right now I don't know.
Bruno Luong
Bruno Luong il 14 Apr 2024 alle 20:12
Ran in:
Another way to find smallest integer x such that A^x = 0 (x is the degree of Nilpotent matrix A):
A =[0 1 2 3 4 5 6 7;
0 0 8 9 10 11 12 13;
0 0 0 1 2 3 4 5;
0 0 0 0 6 7 8 9;
0 0 0 0 0 1 2 3;
0 0 0 0 0 0 4 5;
0 0 0 0 0 0 0 1;
0 0 0 0 0 0 0 0];
J = jordan(A); % symbolictbx required
d0 = diag(J); M main diagonal
isnilpotent = all(d0 == 0);
if isnilpotent
d1 = diag(J,1);
% length of the longest sequence of nonzero above the diagonal
j = diff([0; d1~=0; 0]);
l = diff(find(j));
maxl = max(l);
x = sum([maxl 1]); % it return 1 if dh is empty, add 1 otherwise
else
x = [];
end
x
x = 8

Accedi per commentare.

Più risposte (0)

Categorie

Scopri di più su Linear Algebra in Help Center e File Exchange

Tag

Prodotti


Release

R2024a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by