Zum Inhalt springen

QZ-Algorithmus

aus Wikipedia, der freien Enzyklopädie

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

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

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

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

Vortransformation

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

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

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

Durch analoges spaltenweises Erzeugen von Nullen in <math>A</math> erhält man eine obere Hessenbergmatrix:

  1. <math> A= \begin{pmatrix}*&*&*&*\\{*}&*&*&*\\0&*&*&*\\0&*&*&*\end{pmatrix}, B=\begin{pmatrix}*&*&*&*\\0&*&*&*\\0&*&*&*\\0&0&0&*\end{pmatrix}</math>
  2. <math> A= \begin{pmatrix}*&*&*&*\\{*}&*&*&*\\0&*&*&*\\0&*&*&*\end{pmatrix}, B=\begin{pmatrix}*&*&*&*\\0&*&*&*\\0&0&*&*\\0&0&0&*\end{pmatrix}</math>
  3. <math> A= \begin{pmatrix}*&*&*&*\\{*}&*&*&*\\0&*&*&*\\0&0&*&*\end{pmatrix}, B=\begin{pmatrix}*&*&*&*\\0&*&*&*\\0&0&*&*\\0&0&*&*\end{pmatrix}</math>
  4. <math> A= \begin{pmatrix}*&*&*&*\\{*}&*&*&*\\0&*&*&*\\0&0&*&*\end{pmatrix}, B=\begin{pmatrix}*&*&*&*\\0&*&*&*\\0&0&*&*\\0&0&0&*\end{pmatrix}</math>.

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

QZ-Algorithmus mit impliziten Shifts

1. <math>q:=0</math>

2. while <math> q<n </math> do

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

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

5. Partitioniere <math>B</math> wie <math>A</math>, d. h. <math> B= \begin{pmatrix}B_{11}&B_{12}&B_{13}\\0&B_{22}&B_{23}\\0&0&B_{33}\end{pmatrix} </math>, wobei <math>B_{11}\in\mathbb{R}^{p\times p}, 

 B_{22}\in\mathbb{R}^{m\times m}, B_{33}\in\mathbb{R}^{q\times q} </math> obere Dreiecksmatrizen sind.

6. Bringe <math>A_{33}</math> in obere Schurform: Finde orthogonale <math>Q_{33},Z_{33}</math> so, dass <math>A_{33}:=Q_{33}^{T}A_{33}Z_{33}</math> in Schurform und <math>B_{33}:=Q_{33}^{T}B_{33}Z_{33}</math> obere Dreiecksmatrix ist.

Falls erforderlich: Aufdatieren von <math>Q</math> und <math>Z</math>: <math>Q:=Q\mathrm{diag}(I_p,I_{m},Q_{33})</math>, <math>Z:=Z\mathrm{diag}(I_p,I_{m},Z_{33})</math>.

7. if <math>q<n</math>:

if <math>det(B_{22})=0</math>

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

else

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

end if

8. end if

Wahl der Shifts

9. Bestimme Shifts <math>a, b</math> als Eigenwerte von <math>\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} </math>

10. Bestimme <math> (A_{22}B_{22}^{-1}-aI)(A_{22}B_{22}^{-1}-bI)e_1=\begin{pmatrix}\alpha\\\beta\\\gamma\\0\\\vdots\\0\end{pmatrix}</math>

Der implizite QZ-Schritt

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

Für <math>B_{22}</math> folgt nun: <math>\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} </math>.

Ziel ist es nun, die Dreiecksgestalt von <math>B_{22}</math> durch orthogonale Transformationen (Householder-Spiegelungen) von rechts wiederherzustellen:

12. Finde orthogonales <math>Z_1\in\mathbb{R}^{3\times 3}</math> mit <math>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}</math>. Finde dann orthogonales <math>Z_1'\in\mathbb{R}^{2\times 2}</math>, so dass

<math>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}</math>.

Für <math>A_{22}</math> ergibt sich nun: <math>\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}}</math>. D.h., die Hessenbergstruktur von <math>A_{22}</math> 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 <math>Q_1\in\mathbb{R}^{3\times 3}</math>, <math>Q_1'\in\mathbb{R}^{3\times 3}</math> mit

<math>\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}}</math>. Es werden also nacheinander die Vektoren <math>\begin{pmatrix}a_{21}\\a_{31}\\a_{41}\end{pmatrix}</math> und <math>\begin{pmatrix}a_{32}\\a_{42}\\a_{52}\end{pmatrix}</math>auf <math>\begin{pmatrix}*\\0\\0\end{pmatrix}</math>transformiert.

Die Anwendung der Transformation auf <math>\tilde{B}_{22}</math>von links ergibt jedoch

<math>\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}}</math>, d. h. ein Buckel ist jetzt eine Position tiefer entlang der Diagonalen entstanden.

14. Man wiederhole die Schritte 11–13 so lange, bis <math>A_{22}</math> wieder in oberer Hessenberg- und <math>B_{22}</math> 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 <math>B_{22}</math> an der Diagonalposition j mit Transformationen von links führt zu einem Buckel an der entsprechenden Position in <math>A_{22}</math>. Wird dieser Buckel mit Transformationen von rechts eliminiert, führt das zu einem Buckel in <math>B_{22}</math> an der Diagonalposition j+1 usw.

15. Nach <math>m-2</math> Schritten wird das Ziel erreicht und es ergibt sich <math>Q^{T}_{22}=\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</math>. Analog erhält man

<math>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}')</math>.

Falls <math>(A,B)</math>-invarianten Unterräume benötigt werden, ist es notwendig die Matrizen <math>{Q}</math> und <math>Z</math> aufzudatieren: <math>Q:=Q\mathrm{diag}(I_p,Q_{22},I_q)</math>, <math>Z:=Z\mathrm{diag}(I_p,Z_{22},I_q)</math>

16. end while

Bestimmung der Eigenwerte

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

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.
Abgerufen von „https://wiki-de.moshellshocker.dns64.de/index.php?title=QZ-Algorithmus&oldid=1021613“