QZ-Algorithmus

Der QZ-Algorithmus 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 = λ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 − λB $A,B \in \mathbb{C}^{n\times n}$ Gesucht sind orthogonale Matrizen Q und Z, so dass Q^T(A − λB)Z = T − λ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 − λ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₁₁ von oberer Hessenbergform, A₂₂ von unreduzierter oberer Hessenbergform und A₃₃ 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₃₃ in obere Schurform: Finde orthogonale Q₃₃,Z₃₃ 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: = Qdiag(I_p,I_m,Q₃₃), Z: = Zdiag(I_p,I_m,Z₃₃).

7. if q < n:

if det(B₂₂) = 0

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

else

Führe einen impliziten QZ-Schritt für A₂₂,B₂₂ 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 des unteren 2*2-Blockes von $\begin{pmatrix}a_{m-1,m-2}&a_{m-1,m-1}&a_{m-1,m}\\0&a_{m,m-1}&a_{m,m}\end{pmatrix}\begin{pmatrix}b_{m-2,m-2}&b_{m-2,m-1}&b_{m-2,m}\\0&b_{m-1,m-1}&b_{m-1,m}\\0&0&b_{m,m}\end{pmatrix}^{-1}$

10. Bestimme $(A_{22}-aI)(A_{22}-bI)e_1=\begin{pmatrix}\alpha\\\beta\\\gamma\\0\\\vdots\\0\end{pmatrix}$

Der implizite QZ-Schritt

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

Für B₂₂ 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₂₂ 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}$.

Für A₂₂ ergibt sich nun: ${A_{22}}\mathrm{diag}(Z_1,I_{m-3})={\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}}$.

13. Finde orthogonales $Z_1'\in\mathbb{R}^{2\times 2}$ mit $B_{22}\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₂₂ ergibt sich nun analog zur Multiplikation mit diag(Z₁,I_{m − 3}):

${A_{22}}\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}}$.

14. Man wiederhole die Schritte 11-14 nun für den Vektor $\begin{pmatrix}a_{21}\\a_{31}\\a_{41}\end{pmatrix}$, dann für $\begin{pmatrix}a_{32}\\a_{42}\\a_{52}\end{pmatrix}$ usw. Auf diese Art und Weise "wandern" die Nicht-Null-Elemente unterhalb der Nebendiagonale von A Schritt für Schritt aus der Matrix. Diesen Prozess bezeichnet man auch als "Bulge-Chasing".

15. Nach m − 2 Schritten liegt A wieder in oberer Hessenbergform vor. Für Q₂₂ ergibt sich nun $Q_{22}=\mathrm{diag}(Q_1,I_{m-3})\mathrm{diag}(I_{1},Q_2,I_{m-4}) \cdots \mathrm{diag}(I_{m-3},Q_{m-2})$. Analog erhält man für Z₂₂:

$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: = Qdiag(I_p,Q₂₂,I_q), Z: = Zdiag(I_p,Z₂₂,I_q)

16. end while

Bestimmung der Eigenwerte

In den meisten Fällen konvergiert A im QZ-Algorithmus gegen seine Schur-Form. Es gilt $\Lambda(A,B)=\{\frac{a_{ii}}{b_{ii}}:b_{ii}\ne 0, i\in \{1,\cdots,n\}\}$. Falls ein i existiert, für das a_ii = b_ii = 0 ist, so ist $\Lambda(A,B)=\mathbb{C}$. Falls b_ii = 0 und $a_{ii}\ne 0$, so ist der i-te Eigenwert unendlich.

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 Volume II: Eigensystems. SIAM 2001, ISBN 0-89871-503-2.