QZ-Algorithmus

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

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

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

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

Vortransformation

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

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

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

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

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

QZ-Algorithmus mit impliziten Shifts

1. $q:=0$

2. while $q<n$ do

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

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

5. Partitioniere $B$ wie $A$, d. h. $B={\begin{pmatrix}B_{11}&B_{12}&B_{13}\\0&B_{22}&B_{23}\\0&0&B_{33}\end{pmatrix}}$, wobei $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 $A_{33}$ in obere Schurform: Finde orthogonale $Q_{33},Z_{33}$ so, dass $A_{33}:=Q_{33}^{T}A_{33}Z_{33}$ in Schurform und $B_{33}:=Q_{33}^{T}B_{33}Z_{33}$ obere Dreiecksmatrix ist.

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

7. if $q<n$:

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

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

else

Führe einen impliziten QZ-Schritt für $A_{22},B_{22}$ aus: $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

9. Bestimme Shifts $a,b$ als Eigenwerte von ${\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 $(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

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

Für $B_{22}$ folgt nun: ${\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 $B_{22}$ durch orthogonale Transformationen (Householder-Spiegelungen) von rechts wiederherzustellen:

12. Finde orthogonales $Z_{1}\in \mathbb {R} ^{3\times 3}$ mit $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 $Z_{1}'\in \mathbb {R} ^{2\times 2}$, so dass

$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 $A_{22}$ ergibt sich nun: ${\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 $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 $Q''_{1}\in \mathbb {R} ^{3\times 3}$, $Q_{1}'\in \mathbb {R} ^{3\times 3}$ mit

$\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 ${\begin{pmatrix}a_{21}\\a_{31}\\a_{41}\end{pmatrix}}$ und ${\begin{pmatrix}a_{32}\\a_{42}\\a_{52}\end{pmatrix}}$auf ${\begin{pmatrix}*\\0\\0\end{pmatrix}}$transformiert.

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

$\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 $A_{22}$ wieder in oberer Hessenberg- und $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 $B_{22}$ an der Diagonalposition j mit Transformationen von links führt zu einem Buckel an der entsprechenden Position in $A_{22}$. Wird dieser Buckel mit Transformationen von rechts eliminiert, führt das zu einem Buckel in $B_{22}$ an der Diagonalposition j+1 usw.

15. Nach $m-2$ Schritten wird das Ziel erreicht und es ergibt sich $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

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

16. end while

Bestimmung der Eigenwerte

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

• 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.