matrix - Finding generalized eigenvectors numerically in Matlab -

i have matrix such example (my actual matrices can larger)

a = [-1 -2   -0.5;       0  0.5  0;       0  0   -1]; 

that has 2 linearly-independent eigenvalues (the eigenvalue -1 repeated). obtain complete basis generalized eigenvectors. 1 way know how matlab's jordan function in symbolic math toolbox, i'd prefer designed numeric inputs (indeed, 2 outputs, jordan fails large matrices: "error in mupad command: similarity matrix large."). don't need jordan canonical form, notoriously unstable in numeric contexts, matrix of generalized eigenvectors. there function or combination of functions automates in numerically stable way or must 1 use the generic manual method (how stable such procedure)?

note: "generalized eigenvector," mean non-zero vector can used augment incomplete basis of so-called defective matrix. not mean eigenvectors correspond eigenvalues obtained solving generalized eigenvalue problem using eig or qz (though latter usage quite common, i'd it's best avoided). unless can correct me, don't believe 2 same.

update 1 – 5 months later:

see my answer here how obtain generalized eigenvectors symbolically matrices larger 82-by-82 (the limit test matrix in question).

i'm still interested in numeric schemes (or how such schemes might unstable if they're related calculating jordan form). don't wish blindly implement linear algebra 101 method has been marked duplicate of question it's not numerical algorithm, rather pencil-and paper method used evaluate students (i suppose implemented symbolically however). if can point me either implementation of scheme or numerical analysis of it, i'd interested in that.

update 2 – feb. 2015: of above still true tested in r2014b.

as mentioned in comments, if matrix defective, know eigenvectors/eigenvalue pair want consider identical given tolerance, can proceed example below:

% example matrix a: = [1 0 0 0 0;       3 1 0 0 0;       6 3 2 0 0;       10 6 3 2 0;      15 10 6 3 2] % produce eigenvalues , eigenvectors (not generalized ones) [vecs,vals] = eig(a) 

this should output:

vecs =       0         0         0         0    0.0000      0         0         0    0.2236   -0.2236      0         0    0.0000   -0.6708    0.6708      0    0.0000   -0.0000    0.6708   -0.6708 1.0000   -1.0000    1.0000   -0.2236    0.2236   vals =   2     0     0     0     0  0     2     0     0     0  0     0     2     0     0  0     0     0     1     0  0     0     0     0     1 

where see first 3 eigenvectors identical working precision, 2 last ones. here, must know structure of problem , identify identical eigenvectors of identical eigenvalues. here, eigenvalues identical, know ones consider, , assume corresponding vectors 1-2-3 identical , vectors 4-5. (in practice check norm of differences of eigenvectors , compare tolerance)

now proceed compute generalized eigenvectors, ill-conditioned solve matlab's \, because (a - lambda*i) not full rank. use pseudoinverses:

genvec21 = pinv(a - vals(1,1)*eye(size(a)))*vecs(:,1); genvec22 = pinv(a - vals(1,1)*eye(size(a)))*genvec21; genvec1 = pinv(a - vals(4,4)*eye(size(a)))*vecs(:,4); 

which should give:

genvec21 =     -0.0000     0.0000    -0.0000     0.3333          0  genvec22 =      0.0000    -0.0000     0.1111    -0.2222          0  genvec1 =      0.0745    -0.8832     1.5317     0.6298    -3.5889 

which our other generalized eigenvectors. if check these obtain jordan normal form this:

jordanj = [vecs(:,1) genvec21 genvec22 vecs(:,4) genvec1]; jordanj^-1*a*jordanj 

we obtain:

ans =  2.0000    1.0000    0.0000   -0.0000   -0.0000      0    2.0000    1.0000   -0.0000   -0.0000      0    0.0000    2.0000    0.0000   -0.0000      0    0.0000    0.0000    1.0000    1.0000      0    0.0000    0.0000   -0.0000    1.0000 

which our jordan normal form (with working precision errors).