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

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?

[Bruno Luong]

If I could argue ro me all methods fail. 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.

It is well known that such matrix (Jordan form with non-zero on upper diagonal) are notoriously difficult and cannot be well factorized by EIG and alike.

I wouldn't try to find an explanation logic behind the values in this specific example, eig with 1 or 2 arguments simply use different algorithm, IIRC respectively qr algorithm and algorithm based on Schur decomposition. There is no reason to call eig with eye() as second argument as general rule.

[Torsten]

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]

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

Answer 1

[Bruno Luong]

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.

[Paul]

"Actually I think eig with 2 inputs is more approprite for multiple order eigen value"

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?

Would it matter if A has other, non-repeated eigenvalues?

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)?

Alternatively, is there any reason to ever not use eig(A,I), assuming that accuracy of the result is the only consideration? Can eig(A,I) ever give a less accurate result than eig(A)?

"Torsen have found a good trick here."

Maybe "Torsten's trick" will catch on like "Tony's trick" !

[Bruno Luong]

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]

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