Das Finden generalisierter Eigenvektoren numerisch in Matlab

Question

Ich habe eine Matrix wie in diesem Beispiel (meine tatsächlichen Matrizen können viel größer sein)

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

die nur zwei linear unabhängige Eigenwerte hat (der Eigenwert -1 wird wiederholt). Ich möchte eine komplette Basis mit generalisierten Eigenvektoren erhalten. Eine Möglichkeit, wie ich dies tun kann, ist mit der jordan-Funktion in der Symbolic Math-Toolbox von Matlab, aber ich würde etwas bevorzugen, das für numerische Eingaben konzipiert ist (tatsächlich scheitert jordan bei großen Matrizen mit zwei Ausgängen: "Fehler im MuPAD-Befehl: Ähnlichkeitsmatrix zu groß."). Ich benötige nicht die Jordan-Kanonische Form, die in numerischen Kontexten bekanntermaßen instabil ist, sondern einfach eine Matrix generalisierter Eigenvektoren. Gibt es eine Funktion oder Kombination von Funktionen, die dies auf numerisch stabile Weise automatisiert oder muss man die generische manuelle Methode verwenden (wie stabil ist ein solches Verfahren)?

HINWEIS: Mit "generalisierten Eigenvektor" meine ich einen Nicht-Null-Vektor, der verwendet werden kann, um die unvollständige Basis einer sogenannten defekten Matrix zu ergänzen. Ich meine nicht die Eigenvektoren, die den Eigenwerten entsprechen, die aus der Lösung des generalisierten Eigenwertproblems mit eig oder qz erhalten werden (obwohl diese letztere Verwendung ziemlich verbreitet ist, würde ich sagen, dass es am besten vermieden werden sollte). Es sei denn, jemand kann mich korrigieren, ich glaube nicht, dass die beiden gleich sind.

UPDATE 1 – Fünf Monate später:

Siehe meine Antwort hier um generalisierte Eigenvektoren symbolisch für Matrizen größer als 82x82 zu erhalten (die Grenze für die Testmatrix in dieser Frage).

Ich bin immer noch an numerischen Verfahren interessiert (oder wie solche Verfahren instabil sein könnten, wenn sie alle mit der Berechnung der Jordan-Form zusammenhängen). Ich möchte nicht blind das lineare Algebra 101-Verfahren implementieren, das als Duplikat dieser Frage markiert wurde, da es kein numerischer Algorithmus ist, sondern vielmehr eine Bleistift-und-Papier-Methode, die zur Bewertung von Studenten verwendet wird (ich nehme an, dass sie jedoch symbolisch implementiert werden könnte). Wenn jemand mir entweder eine Implementierung dieses Schemas zeigen kann oder eine numerische Analyse davon, wäre ich daran interessiert.

UPDATE 2 – Feb. 2015: Alles oben Genannte ist immer noch gültig, wie in R2014b getestet.

Answer 1

Wie in meinen Kommentaren erwähnt, wenn Ihre Matrix fehlerhaft ist, Sie jedoch wissen, welches Eigenvektor/Eigenwert-Paar Sie unter Berücksichtigung Ihrer Toleranz als identisch betrachten möchten, können Sie wie in diesem Beispiel unten vorgehen:

% Beispiel Matrix A:
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]
% Erzeugen von Eigenwerten und Eigenvektoren (nicht generalisierte)
[vecs,vals] = eig(A)

Dies sollte ausgegeben werden:

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

Hier sehen wir, dass die ersten drei Eigenvektoren fast identisch mit Arbeitsgenauigkeit sind, ebenso wie die beiden letzten. Hier müssen Sie die Struktur Ihres Problems kennen und die identischen Eigenvektoren identischer Eigenwerte identifizieren. Hier sind die Eigenwerte genau identisch, sodass wir wissen, welche betrachtet werden sollen, und wir nehmen an, dass die entsprechenden Vektoren 1-2-3 identisch sind und Vektoren 4-5. (In der Praxis überprüfen Sie wahrscheinlich die Norm der Differenzen der Eigenvektoren und vergleichen sie mit Ihrer Toleranz)

Jetzt gehen wir dazu über, die generalisierten Eigenvektoren zu berechnen, aber dies ist schlecht bedingt, um sie einfach mit dem \ von Matlab zu lösen, denn offensichtlich ist (A - lambda*I) nicht voller Rang. Daher verwenden wir Pseudoinversen:

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

Was sollte geben:

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

Dies sind unsere anderen generalisierten Eigenvektoren. Wenn wir sie jetzt überprüfen, um die Jordan-Normalform wie folgt zu erhalten:

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

Erhalten wir:

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

Was unsere Jordan-Normalform ist (mit Arbeitsgenauigkeitsfehlern).