QZ-Algorithmus

Der QZ-Algorithmus oder die QZ-Faktorisierung ist ein numerisches Verfahren zur Lösung des verallgemeinerten Eigenwertproblems.

${\displaystyle Ax=\lambda Bx^{\,}}$ , mit ${\displaystyle A,B\in \mathbb {R} ^{n\times n}}$ bzw. ${\displaystyle A,B\in \mathbb {C} ^{n\times n}}$

Das verallgemeinerte Eigenwertproblem ist äquivalent zum Eigenwertproblem ${\displaystyle AB^{-1}y=\lambda y}$, wobei ${\displaystyle y=Bx}$ und ${\displaystyle B}$ invertierbar sein muss. Es wird jedoch nicht explizit die Matrix ${\displaystyle B^{-1}}$ berechnet, um die Kondition des Problems nicht zu verschlechtern, sondern ${\displaystyle A}$ und ${\displaystyle B}$ werden simultan durch Ähnlichkeitstransformationen (Givens-Rotationen und Householder-Spiegelungen) in verallgemeinerte Schurform gebracht.

Gegeben ist ein Matrixbüschel ${\displaystyle A-\lambda B}$. Gesucht sind orthogonale Matrizen ${\displaystyle Q}$ und ${\displaystyle Z}$, so dass ${\displaystyle Q^{T}(A-\lambda B)Z=T-\lambda S}$ von verallgemeinerter Schurform ist, d. h. ${\displaystyle T}$ ist von quasi-oberer Dreiecksform und ${\displaystyle S}$ ist von oberer Dreiecksform. Im Fall ${\displaystyle A,B\in \mathbb {C} ^{n\times n}}$ ist ${\displaystyle T}$ stets von oberer Dreiecksform. Aus der verallgemeinerten Schurform lassen sich dann die Eigenwerte und aus ${\displaystyle Q}$ und ${\displaystyle Z}$ ${\displaystyle (A,B)}$-invariante Unterräume des Matrixbüschels ${\displaystyle A-\lambda B}$ bestimmen.

Vortransformation[Bearbeiten | Quelltext bearbeiten]

Ziel dieses Schrittes ist es, die Matrix ${\displaystyle A}$ durch orthogonale Transformationen auf obere Hessenbergform und die Matrix ${\displaystyle B}$ auf obere Dreiecksform zu bringen. Durch ${\displaystyle n-1}$ Householder-Spiegelungen von links wird ${\displaystyle B}$ auf obere Dreiecksform transformiert. Wendet man die gleichen Transformationen gleichzeitig auf ${\displaystyle A}$ an, ergibt sich (Veranschaulichung an einem Beispiel der Größe (4,4)): ${\displaystyle A={\begin{pmatrix}*&*&*&*\\{*}&*&*&*\\{*}&*&*&*\\{*}&*&*&*\end{pmatrix}},B={\begin{pmatrix}*&*&*&*\\0&*&*&*\\0&0&*&*\\0&0&0&*\end{pmatrix}}}$.

Man finde nun eine Givens-Rotation, die von links angewendet auf A folgende Matrix ergibt: ${\displaystyle A={\begin{pmatrix}*&*&*&*\\{*}&*&*&*\\{*}&*&*&*\\0&*&*&*\end{pmatrix}}}$. Damit erhält man für ${\displaystyle B={\begin{pmatrix}*&*&*&*\\0&*&*&*\\0&0&*&*\\0&0&*&*\end{pmatrix}}}$.

Durch Anwendung einer Givens-Rotation von rechts kann die obere Dreiecksform von ${\displaystyle B}$ wiederhergestellt werden, ohne die Null an der linken unteren Position von A zu zerstören: ${\displaystyle A={\begin{pmatrix}*&*&*&*\\{*}&*&*&*\\{*}&*&*&*\\0&*&*&*\end{pmatrix}},B={\begin{pmatrix}*&*&*&*\\0&*&*&*\\0&0&*&*\\0&0&0&*\end{pmatrix}}}$.

Durch analoges spaltenweises Erzeugen von Nullen in ${\displaystyle A}$ erhält man eine obere Hessenbergmatrix:

1. ${\displaystyle A={\begin{pmatrix}*&*&*&*\\{*}&*&*&*\\0&*&*&*\\0&*&*&*\end{pmatrix}},B={\begin{pmatrix}*&*&*&*\\0&*&*&*\\0&*&*&*\\0&0&0&*\end{pmatrix}}}$
2. ${\displaystyle A={\begin{pmatrix}*&*&*&*\\{*}&*&*&*\\0&*&*&*\\0&*&*&*\end{pmatrix}},B={\begin{pmatrix}*&*&*&*\\0&*&*&*\\0&0&*&*\\0&0&0&*\end{pmatrix}}}$
3. ${\displaystyle A={\begin{pmatrix}*&*&*&*\\{*}&*&*&*\\0&*&*&*\\0&0&*&*\end{pmatrix}},B={\begin{pmatrix}*&*&*&*\\0&*&*&*\\0&0&*&*\\0&0&*&*\end{pmatrix}}}$
4. ${\displaystyle A={\begin{pmatrix}*&*&*&*\\{*}&*&*&*\\0&*&*&*\\0&0&*&*\end{pmatrix}},B={\begin{pmatrix}*&*&*&*\\0&*&*&*\\0&0&*&*\\0&0&0&*\end{pmatrix}}}$.

Falls ${\displaystyle (A,B)}$-invariante Unterräume berechnet werden sollen, so ist es notwendig, das Produkt der Transformationsmatrizen, die jeweils von links auf ${\displaystyle A}$ und ${\displaystyle B}$ angewendet werden, in einer Matrix ${\displaystyle Q}$ und das Produkt der Transformationsmatrizen, die von rechts angewendet werden, in einer Matrix ${\displaystyle Z}$ zu speichern.

QZ-Algorithmus mit impliziten Shifts[Bearbeiten | Quelltext bearbeiten]

1. ${\displaystyle q:=0}$

2. while ${\displaystyle q<n}$ do

3. Bestimme alle ${\displaystyle j\in \{1,\cdots ,n-1\}}$ mit ${\displaystyle |a_{j+1,j}|\leq \varepsilon (|a_{j,j}|+|a_{j+1,j+1}|)}$. Für diese ${\displaystyle j}$ setze ${\displaystyle a_{j,j+1}=0}$.

4. Deflation: Finde minimales ${\displaystyle p}$ und maximales ${\displaystyle q}$ mit ${\displaystyle p,q\in \{1,\cdots ,n\}}$ und definiere ${\displaystyle m:=n-p-q}$, so dass gilt: ${\displaystyle A={\begin{pmatrix}A_{11}&A_{12}&A_{13}\\0&A_{22}&A_{23}\\0&0&A_{33}\end{pmatrix}}}$, wobei ${\displaystyle A_{11}\in \mathbb {R} ^{p\times p},A_{22}\in \mathbb {R} ^{m\times m},A_{33}\in \mathbb {R} ^{q\times q}}$ und ${\displaystyle A_{11}}$ von oberer Hessenbergform, ${\displaystyle A_{22}}$ von unreduzierter oberer Hessenbergform und ${\displaystyle A_{33}}$ von quasi-oberer Dreiecksform ist.

5. Partitioniere ${\displaystyle B}$ wie ${\displaystyle A}$, d. h. ${\displaystyle B={\begin{pmatrix}B_{11}&B_{12}&B_{13}\\0&B_{22}&B_{23}\\0&0&B_{33}\end{pmatrix}}}$, wobei ${\displaystyle B_{11}\in \mathbb {R} ^{p\times p},B_{22}\in \mathbb {R} ^{m\times m},B_{33}\in \mathbb {R} ^{q\times q}}$ obere Dreiecksmatrizen sind.

6. Bringe ${\displaystyle A_{33}}$ in obere Schurform: Finde orthogonale ${\displaystyle Q_{33},Z_{33}}$ so, dass ${\displaystyle A_{33}:=Q_{33}^{T}A_{33}Z_{33}}$ in Schurform und ${\displaystyle B_{33}:=Q_{33}^{T}B_{33}Z_{33}}$ obere Dreiecksmatrix ist.

Falls erforderlich: Aufdatieren von ${\displaystyle Q}$ und ${\displaystyle Z}$: ${\displaystyle Q:=Q\mathrm {diag} (I_{p},I_{m},Q_{33})}$, ${\displaystyle Z:=Z\mathrm {diag} (I_{p},I_{m},Z_{33})}$.

7. if ${\displaystyle q<n}$:

if ${\displaystyle det(B_{22})=0}$

Transformiere mithilfe einer Givens-Rotation von rechts ${\displaystyle a_{n-q,n-q-1}=0}$, um die Rang-Defizienz von ${\displaystyle B_{22}}$ auf ${\displaystyle B_{33}}$ zu verschieben. Durch die Annullierung von ${\displaystyle a_{n-q,n-q-1}}$ ist ${\displaystyle A_{22}}$ keine unreduzierte Hessenbergmatrix mehr, somit wird ${\displaystyle q}$ erhöht und es besteht die Möglichkeit, dass ${\displaystyle B_{22}}$ in der neuen Partionierung regulär ist.

else

Führe einen impliziten QZ-Schritt für ${\displaystyle A_{22},B_{22}}$ aus: ${\displaystyle A_{22}:=Q_{22}^{T}A_{22}Z_{22},\quad {B_{22}}:=Q_{22}^{T}{B_{22}}Z_{22}}$.

end if

8. end if

Wahl der Shifts[Bearbeiten | Quelltext bearbeiten]

9. Bestimme Shifts ${\displaystyle a,b}$ als Eigenwerte von ${\displaystyle {\begin{pmatrix}a_{m-1,m-1}&a_{m-1,m}\\a_{m,m-1}&a_{m,m}\end{pmatrix}}{\begin{pmatrix}b_{m-1,m-1}&b_{m-1,m}\\0&b_{m,m}\end{pmatrix}}^{-1}}$

10. Bestimme ${\displaystyle (A_{22}B_{22}^{-1}-aI)(A_{22}B_{22}^{-1}-bI)e_{1}={\begin{pmatrix}\alpha \\\beta \\\gamma \\0\\\vdots \\0\end{pmatrix}}}$

Der implizite QZ-Schritt[Bearbeiten | Quelltext bearbeiten]

11. Finde orthogonales ${\displaystyle Q_{1}}$ mit ${\displaystyle Q_{1}^{T}{\begin{pmatrix}\alpha \\\beta \\\gamma \end{pmatrix}}={\begin{pmatrix}*\\0\\0\end{pmatrix}}}$

Für ${\displaystyle B_{22}}$ folgt nun: ${\displaystyle {\begin{pmatrix}Q_{1}^{T}&0\\0&I_{m-3}\end{pmatrix}}B_{22}={\begin{pmatrix}*&*&*&\cdots &\cdots &*\\{*}&*&*&\cdots &\cdots &*\\{*}&*&*&\cdots &\cdots &*\\0&0&0&\ddots &&\vdots \\\vdots &\vdots &\vdots &&\ddots &\vdots \\0&0&0&\cdots &\cdots &*\end{pmatrix}}}$.

Ziel ist es nun, die Dreiecksgestalt von ${\displaystyle B_{22}}$ durch orthogonale Transformationen (Householder-Spiegelungen) von rechts wiederherzustellen:

12. Finde orthogonales ${\displaystyle Z_{1}\in \mathbb {R} ^{3\times 3}}$ mit ${\displaystyle B_{22}\mathrm {diag} (Z_{1},I_{m-3})={\begin{pmatrix}*&*&*&\cdots &\cdots &*\\{*}&*&*&\cdots &\cdots &*\\0&0&*&\cdots &\cdots &*\\0&0&0&\ddots &&\vdots \\\vdots &\vdots &\vdots &&\ddots &\vdots \\0&0&0&\cdots &\cdots &*\end{pmatrix}}}$. Finde dann orthogonales ${\displaystyle Z_{1}'\in \mathbb {R} ^{2\times 2}}$, so dass

${\displaystyle B_{22}\mathrm {diag} (Z_{1},I_{m-3})\mathrm {diag} (Z_{1}',I_{m-2})={\begin{pmatrix}*&*&*&\cdots &\cdots &*\\0&*&*&\cdots &\cdots &*\\0&0&*&\cdots &\cdots &*\\0&0&0&\ddots &&\vdots \\\vdots &\vdots &\vdots &&\ddots &\vdots \\0&0&0&\cdots &\cdots &*\end{pmatrix}}}$.

Für ${\displaystyle A_{22}}$ ergibt sich nun: ${\displaystyle {\tilde {A}}_{22}:={A_{22}}\mathrm {diag} (Z_{1},I_{m-3})\mathrm {diag} (Z'_{1},I_{m-2})={\begin{pmatrix}*&*&*&\cdots &\cdots &\cdots &*\\{*}&*&*&\cdots &\cdots &\cdots &*\\{*}&*&*&\cdots &\cdots &\cdots &*\\{*}&*&*&\cdots &\cdots &\cdots &*\\0&0&0&\ddots &&&\vdots \\\vdots &\vdots &\vdots &&\ddots &&\vdots \\0&0&0&\cdots &0&*&*\end{pmatrix}}}$. D.h., die Hessenbergstruktur von ${\displaystyle A_{22}}$ ist durch einen unerwünschten 2x2 "Buckel" zerstört.

13. Dieser Buckel kann durch elementäre, orthogonale Transformationen (z. B. Householder-Spiegelungen) von links eliminiert werden. Finde also orthogonales ${\displaystyle Q''_{1}\in \mathbb {R} ^{3\times 3}}$, ${\displaystyle Q_{1}'\in \mathbb {R} ^{3\times 3}}$ mit

${\displaystyle \mathrm {diag} (1,Q'_{1},I_{m-4})^{T}\mathrm {diag} (I_{2},Q''_{1},I_{m-5})^{T}{{\tilde {A}}_{22}}={\begin{pmatrix}*&*&*&\cdots &\cdots &\cdots &*\\*&*&*&\cdots &\cdots &\cdots &*\\0&*&*&\ddots &\cdots &\cdots &*\\0&0&*&\ddots &\cdots &\cdots &*\\0&0&0&*&&\vdots \\\vdots &\vdots &\vdots &&\ddots &&\vdots \\0&0&0&\cdots &0&*&*\end{pmatrix}}}$. Es werden also nacheinander die Vektoren ${\displaystyle {\begin{pmatrix}a_{21}\\a_{31}\\a_{41}\end{pmatrix}}}$ und ${\displaystyle {\begin{pmatrix}a_{32}\\a_{42}\\a_{52}\end{pmatrix}}}$auf ${\displaystyle {\begin{pmatrix}*\\0\\0\end{pmatrix}}}$transformiert.

Die Anwendung der Transformation auf ${\displaystyle {\tilde {B}}_{22}}$von links ergibt jedoch

${\displaystyle \mathrm {diag} (1,Q'_{1},I_{m-4})^{T}\mathrm {diag} (I_{2},Q''_{1},I_{m-5})^{T}{{\tilde {B}}_{22}}={\begin{pmatrix}*&*&*&\cdots &\cdots &\cdots &*\\0&*&*&\cdots &\cdots &\cdots &*\\0&*&*&\ddots &\cdots &\cdots &*\\0&*&*&*&\ddots &\cdots &*\\0&0&0&0&*&\vdots \\\vdots &\vdots &\vdots &&\ddots &&\vdots \\0&0&0&\cdots &0&0&*\end{pmatrix}}}$, d. h. ein Buckel ist jetzt eine Position tiefer entlang der Diagonalen entstanden.

14. Man wiederhole die Schritte 11–13 so lange, bis ${\displaystyle A_{22}}$ wieder in oberer Hessenberg- und ${\displaystyle B_{22}}$ wieder in oberer Dreieckstruktur vorliegt. Diesen Prozess bezeichnet man, analog zum QR-Algorithmus, auch als "Buckel-Jagen" oder "Bulge-Chasing". Die Eliminierung eines Buckels in ${\displaystyle B_{22}}$ an der Diagonalposition j mit Transformationen von links führt zu einem Buckel an der entsprechenden Position in ${\displaystyle A_{22}}$. Wird dieser Buckel mit Transformationen von rechts eliminiert, führt das zu einem Buckel in ${\displaystyle B_{22}}$ an der Diagonalposition j+1 usw.

15. Nach ${\displaystyle m-2}$ Schritten wird das Ziel erreicht und es ergibt sich ${\displaystyle Q_{22}^{T}=\mathrm {diag} (Q_{1},I_{m-3})^{T}\mathrm {diag} (1,Q'_{1},I_{m-4})^{T}\mathrm {diag} (I_{2},Q_{1}'',I_{m-5})^{T}\cdots \mathrm {diag} (I_{m-3},Q_{m-2})^{T}}$. Analog erhält man

${\displaystyle Z_{22}=\mathrm {diag} (Z_{1},I_{m-3})\mathrm {diag} (Z_{1}',I_{m-2})\cdots \mathrm {diag} (I_{m-2},Z_{m-2})\mathrm {diag} (I_{m-2},Z_{m-2}')}$.

Falls ${\displaystyle (A,B)}$-invarianten Unterräume benötigt werden, ist es notwendig die Matrizen ${\displaystyle {Q}}$ und ${\displaystyle Z}$ aufzudatieren: ${\displaystyle Q:=Q\mathrm {diag} (I_{p},Q_{22},I_{q})}$, ${\displaystyle Z:=Z\mathrm {diag} (I_{p},Z_{22},I_{q})}$

16. end while

Bestimmung der Eigenwerte[Bearbeiten | Quelltext bearbeiten]

In den meisten Fällen konvergiert ${\displaystyle (A,B)}$ im QZ-Algorithmus gegen seine verallgemeinerte, reelle Schur-Form. Für skalare Diagonalblöcke in A gilt ${\displaystyle \lambda _{i}={\frac {a_{ii}}{b_{ii}}}:b_{ii}\neq 0}$ und ${\displaystyle \lambda _{i}=\infty }$ falls ${\displaystyle b_{ii}=0}$. Falls ein ${\displaystyle i}$ existiert, für das ${\displaystyle a_{ii}=b_{ii}=0}$ ist, so ist ${\displaystyle \Lambda (A,B)=\mathbb {C} }$. ${\displaystyle 2\times 2}$ Diagonalblöcke von ${\displaystyle A}$ beziehen sich (analog zum QR-Algorithmus) auf Paare komplex konjugierter Eigenwerte ${\displaystyle \lambda ,{\overline {\lambda }}=\Lambda \left({\begin{pmatrix}a_{ii}&a_{i,i+1}\\a_{i+1,i}&a_{i+1,i+1}\end{pmatrix}},{\begin{pmatrix}b_{ii}&b_{i,i+1}\\0&b_{i+1,i+1}\end{pmatrix}}\right)}$.

Literatur[Bearbeiten | Quelltext bearbeiten]

• Gene H. Golub, Charles F. Van Loan: Matrix Computations. Johns Hopkins University Press, 1996, ISBN 0-8018-5414-8.
• G. W. Stewart: Matrix Algorithms. Band II: Eigensystems. SIAM 2001, ISBN 0-89871-503-2.