QR-Algorithmus
Der QR-Algorithmus ist ein numerisches Verfahren zur Berechnung aller Eigenwerte und eventuell der Eigenvektoren einer quadratischen Matrix. Das auch QR-Verfahren oder QR-Iteration genannte Verfahren basiert auf der QR-Zerlegung und wurde im Jahre 1961–1962 unabhängig voneinander von John G. F. Francis und Wera Nikolajewna Kublanowskaja eingeführt. Ein Vorläufer war der LR-Algorithmus von Heinz Rutishauser (1958), der aber weniger stabil ist und auf der LR-Zerlegung basiert. Oft konvergieren die Iterierten aus dem QR-Algorithmus gegen die Schur-Form der Matrix. Das originale Verfahren ist recht aufwendig und damit – selbst auf heutigen Rechnern – für Matrizen mit hunderttausenden Zeilen und Spalten nicht praktikabel.
Abgeleitete Varianten wie das Multishift-Verfahren von Z. Bai und James Demmel 1989 und die numerisch stabilere Variante von K. Braman, R. Byers und R. Mathias 2002 haben praktische Laufzeiten, die kubisch in der Größe der Matrix sind. Letzteres Verfahren ist in der numerischen Softwarebibliothek LAPACK implementiert, die wiederum in vielen Computeralgebrasystemen (CAS) für die numerischen Matrixalgorithmen verwendet wird.
Beschreibung des QR-Algorithmus
Ziel des Rechenverfahrens
Ausgehend von einer reellen oder komplexen Matrix
bzw.
wird eine Folge orthogonaler oder unitärer Matrizen
bestimmt, so dass die rekursive Matrixfolge
weitestgehend gegen eine obere Dreiecksmatrix konvergiert. Da alle Umformungen
in der Rekursion Ähnlichkeitstransformationen sind, haben alle Matrizen der
Matrixfolge
dieselben Eigenwerte mit denselben Vielfachheiten.
Wird im Grenzwert eine obere Dreiecksmatrix erreicht, eine Schurform
von ,
so lassen sich die Eigenwerte aus den Diagonalelementen ablesen. Im Falle einer
reellen Matrix mit orthogonalen Umformungen kann es jedoch auch komplexe
Eigenwerte geben. Dann ist das Ergebnis des Verfahrens eine obere
Blockdreiecksmatrix, deren Diagonalblöcke die Größe
für die reellen Eigenwerte oder die Größe
für komplexe Eigenwerte haben. Einem Eigenwert
und seinem konjugiert komplexen Eigenwert entspricht dabei der Diagonalblock der
entsprechenden Drehstreckung
.
Allgemeines Schema des Verfahrens
Ausgehend von einer Matrix
(bzw.
)
wird eine Folge von Matrizen
nach folgender Vorschrift bestimmt:
- for
do
- Bestimme ein Polynom
und werte es in der Matrix
aus.
- Berechne die QR-Zerlegung
von
.
- Berechne die neue Iterierte:
.
- end for
In dieser allgemeinen Form können auch die Varianten mit (impliziten) Shifts (engl. für Verschiebung) und Mehrfach-Shifts dargestellt und analysiert werden.
Die Matrixfunktion
ist oft ein Polynom (hier also eine Linearkombination von Matrixpotenzen) mit
reellen (bzw. komplexen) Koeffizienten. Im einfachsten Fall wird ein lineares
Polynom gewählt, das dann die Gestalt
hat. Der Algorithmus vereinfacht sich in diesem Fall auf die klassische
Version (mit Einfach-Shift):
- for
do
- Bestimme eine geeignete Zahl
- Berechne die QR-Zerlegung
von
- Berechne die neue Iterierte:
- end for
Wird in jedem Iterationsschritt
gesetzt, so stimmt das Verfahren mit der auf Unterräume (hier den vollen
Vektorraum) erweiterten Potenzmethode
überein.
Das die QR-Zerlegung steuernde Polynom kann von Anfang an fixiert sein oder
dynamisch in jedem Iterationsschritt der aktuellen Matrix
angepasst werden. Es gibt auch Versuche, für
rationale Matrixfunktionen zu verwenden, d.h. Funktionen der Form
mit Polynomen
und
.
Deflation
Ergibt es sich in einem Iterationsschritt, dass ein linksunterer Block aus
den Spalten
und deren Zeilen
in den Beträgen aller seiner Komponententen eine vorher bestimmte Schwelle
unterschreitet, so wird das Verfahren auf den zwei diagonalen quadratischen
Teilblöcken der Zeilen/Spalten
sowie
separat fortgesetzt. Sind die durch Teilung entstandenen Matrizen klein genug,
so kann die Bestimmung der Eigenwerte z.B. mittels Berechnung des
charakteristischen Polynoms und dessen Nullstellen beendet werden.
Transformation auf Hessenberg-Form
Das Ziel des QR-Algorithmus ist es, eine obere Dreiecksform, oder eine
Block-Dreiecksform herzustellen, in der Blöcke der Größe
komplexen Eigenwerten entsprechen. Das kann nahezu auf einfache Weise durch eine
Ähnlichkeitstransformation auf Hessenberg-Form,
d.;h. auf eine Matrix mit nur einer (nicht identisch verschwindenden)
unteren Nebendiagonalen, erreicht werden.
- Setze
- für k=1 bis n-1 führe aus
-
- Bilde das Spaltensegment
- Bestimme die Householder-Spiegelung
, die
auf den ersten kanonischen Basisvektor abbildet,
- Führe mit der Blockmatrix
die Ersetzung von
durch die ähnliche Matrix
durch.
- Bilde das Spaltensegment
- Vermerke die Gesamttransformationsmatrix
.
befindet sich nun in Hessenberg-Form.
Durch die Hessenberg-Form wird die Bestimmung der charakteristischen Polynome von Teilmatrizen erleichtert. Die Hessenberg-Form einer symmetrischen Matrix ist eine Tridiagonalform, was weitere Rechnungen wesentlich beschleunigt.
Weiterhin wird in jedem Schritt des QR-Algorithmus die Hessenberg-Form erhalten. Grundlage hierfür ist die Arithmetik verallgemeinerter Dreiecksmatrizen, bei denen alle Einträge unterhalb einer unteren Nebendiagonalen verschwinden. Eine Hessenberg-Matrix ist demnach eine verallgemeinerte Dreiecksmatrix mit einer Nebendiagonalen. Unter Multiplikation addieren sich die Anzahlen nichtverschwindender Nebendiagonalen, bei Addition bleibt meist die größere Anzahl erhalten.
Daher haben
wie auch die orthogonale Matrix
die Anzahl von
unteren Nebendiagonalen. Nun gilt, wegen
,
auch
,
und letzteres Produkt hat ebenfalls
Nebendiagonalen. Das ist im Allgemeinen nur möglich, wenn
genau eine Nebendiagonale aufweist, also wieder in Hessenbergform ist. Aus der
Zerlegung von
in Linearfaktoren folgt (s. unten), dass dies immer der Fall ist.
Man kann darauf aufbauend zeigen (Francis 1962 nach Bai/Demmel), dass schon
die erste Spalte
von
,
die auch proportional zur ersten Spalte von
ist, die nachfolgende Matrix
vollständig bestimmt. Genauer, ist
eine orthogonale Matrix, deren erste Spalte proportional zu
ist, so entsteht
,
indem die transformierte Matrix
wieder in Hessenbergform gebracht wird. Da in
nur die Komponenten der Zeilen
von Null verschieden sind, kann
als eine Modifikation der Einheitsmatrix
im linksoberen
-Block
sein, mit einem
.
Varianten des QR-Algorithmus
Einfache QR-Iteration
Es wird
gewählt. Das Verfahren kann dann kurz als QR-Zerlegung
gefolgt vom Zusammenmultiplizieren
in umgekehrter Reihenfolge angegeben werden. Dieses Verfahren ist die direkte
Verallgemeinerung der simultanen Potenzmethode zur Bestimmung der ersten
Eigenwerte einer Matrix. Dieser Zusammenhang wird bei der Unterraumiteration
hergeleitet. Indirekt wird auch eine simultane inverse Potenzmethode
ausgeführt.
QR-Algorithmus mit einfachen Shifts
Es wird
gesetzt. Damit kann das Verfahren alternativ auch als QR-Zerlegung
gefolgt vom um den Shift korrigierten Zusammenmultiplizieren
dargestellt werden. Üblicherweise wird versucht, mit dem Shift
den betragskleinsten Eigenwert zu approximieren. Dazu kann das letzte
Diagonalelement
gewählt werden. Die einfache QR-Iteration ergibt sich, indem alle Shifts zu Null
gesetzt werden.
Besitzt
Hessenberg-Form, so muss
als Produkt
einer Matrix mit und einer Matrix ohne Nebendiagonalen
ebenfalls Hessenberg-Form besitzen. Desgleichen gilt daher auch für
.
Wird also in Vorbereitung des QR-Algorithmus in
auf Hessenberg-Form gebracht, so bleibt dies während des gesamten Algorithmus
erhalten.
Einfache Shifts für symmetrische Matrizen
Eine symmetrische reelle Matrix
hat ausschließlich reelle Eigenwerte. Die Symmetrie bleibt während des
QR-Algorithmus in allen
erhalten. Für symmetrische Matrizen schlug Wilkinson (1965) vor, denjenigen
Eigenwert der rechtsuntersten
-Teilmatrix
als Shift zu wählen, der näher an
liegt. Wilkinson zeigte, dass die so bestimmte Matrixfolge
gegen eine Diagonalmatrix konvergiert, deren Diagonalelemente die Eigenwerte von
sind. Die Konvergenzgeschwindigkeit ist dabei quadratisch.
QR-Algorithmus mit doppelten Shifts
Es kann ein Paar von einfachen Shifts in einem Iterationsschritt zusammengefasst werden. In der Konsequenz bedeutet dies, dass für reelle Matrizen auf komplexe Shifts verzichtet werden kann. In der oben angegebenen Notation ist
eine QR-Zerlegung für das quadratische Polynom ,
ausgewertet in
.
Die Koeffizienten dieses Polynoms sind auch für ein konjugiertes Paar komplexer
Shifts reell. Somit können auch komplexe Eigenwerte reeller Matrizen
approximiert werden, ohne dass in der Rechnung komplexe Zahlen auftauchen.
Eine übliche Wahl für diesen Doppelshift besteht aus den Eigenwerten der
rechtsuntersten -Teilmatrix,
d.h. das quadratische Polynom ist das charakteristische dieses Blocks,
.
QR-Algorithmus mit multiplen Shifts
Es wird eine Zahl
größer
,
aber wesentlich kleiner als die Größe
der Matrix
festgelegt. Das Polynom
kann als das charakteristische Polynom der rechtsuntersten quadratischen
-Teilmatrix
der aktuellen Matrix
gewählt werden. Eine andere Strategie besteht darin, die
Eigenwerte der rechtsuntersten
-Teilmatrix
zu bestimmen und die
betragskleinsten Eigenwerte
darunter auszuwählen. Mit diesen wird dann eine QR-Zerlegung von
und
bestimmt.
Mit einem Multishift-Verfahren wird oft erreicht, dass die Komponenten des
linksunteren -Blocks
in der Folge der iterierten Matrizen besonders schnell klein werden und somit
eine Aufspaltung des Eigenwertproblems erreicht wird.
Implizite Multishift-Iteration
Das Zusammenfassen mehrfacher Shifts in der allgemeinen Form ist sehr
aufwendig. Wie oben angesprochen, kann der Aufwand dadurch verringert werden,
dass in einem vorbereitenden Schritt in
die Hessenberg-Form hergestellt wird. Da jeder Multishift-Schritt aus einzelnen
(auch komplexen) Shifts zusammengesetzt werden kann, bleibt die Hessenberg-Form
während des gesamten Algorithmus erhalten.
Dadurch kann der QR-Algorithmus in einen „Bulge-chasing“-Algorithmus umgewandelt werden, der eine Delle in der Hessenbergform am oberen Diagonalenende erzeugt und diese dann die Diagonale herunter und am unteren Ende aus der Matrix „jagt“.
- for
do
- Berechne das Polynom
nach einer der angegebenen Varianten,
- Bestimme den Vektor
.
- Bestimme eine Spiegelung von
auf den ersten Einheitsvektor. Da in
nur die ersten
Komponenten nicht verschwinden, hat diese Spiegelung eine einfache Blockgestalt.
- Bilde die Matrix
und transformiere sie so, dass
wieder in Hessenberg-Form ist.
- end for
Wird
aus Householder-Spiegelungen zusammengesetzt,
,
so haben diese Blockdiagonalgestalt
.
Anmerkungen zur Funktionsweise
Ähnlichkeitstransformationen
Die im QR-Algorithmus berechneten Matrizen sind zueinander unitär ähnlich, da
aufgrund von
gilt. Damit haben alle Matrizen
dieselben Eigenwerte (mit der algebraischen und geometrischen Vielfachheit
gezählt).
Wahl der Shifts, Konvergenz
Eine einfache aber nicht sehr effektive Wahl ist die Wahl von Shifts
identisch Null. Die Iterierten
des resultierenden Algorithmus, des QR-Algorithmus in der Grundform,
konvergieren im Wesentlichen, wenn alle Eigenwerte dem Betrage nach verschieden
sind, gegen eine obere Dreiecksmatrix mit den Eigenwerten auf der Diagonalen.
Konvergenz im Wesentlichen bedeutet, dass die Elemente des unteren Dreiecks von
gegen Null gehen und die Diagonalelemente gegen die Eigenwerte. Über die
Elemente im oberen Dreieck wird also nichts ausgesagt.
Werden die Shifts als das Matrixelement unten rechts der aktuellen Iterierten
gewählt, also ,
so konvergiert der Algorithmus unter geeigneten Umständen quadratisch oder sogar
kubisch. Dieser Shift ist als Rayleigh-Quotienten-Shift bekannt, da er über die
Inverse
Iteration mit einem Rayleigh-Quotienten
im Zusammenhang steht.
Bei der Rechnung im Reellen ()
möchte man die reelle Schur-Form berechnen und trotzdem mit konjugiert komplexen
Eigenwerten arbeiten können. Dazu gibt es verschiedene Shift-Strategien.
Eine Erweiterung von einfachen Shifts ist der nach James Hardy
Wilkinson benannte Wilkinson-Shift. Für den Wilkinson-Shift wird der näher
am letzten Matrixelement liegende Eigenwert der letzten
Hauptunterabschnittsmatrix
verwendet.
Der QR-Algorithmus als Erweiterung der Potenzmethode
Zur Analyse des Algorithmus werden die zusätzlichen Matrixfolgen der
kumulierten Produkte
und
,
definiert. Dabei sind die Produkte
von orthogonalen bzw. unitären Matrizen wieder orthogonale bzw. unitäre
Matrizen, genauso sind die Produkte
von rechtsoberen Dreiecksmatrizen wieder rechtsobere Dreiecksmatrizen. Die
Matrizen der QR-Iteration ergeben sich durch Ähnlichkeitstransformationen aus
,
denn
.
Daraus folgert man auf QR-Zerlegungen der Potenzen von :
Es werden also im Verlaufe des Algorithmus implizit QR-Zerlegungen der
Potenzen der Matrix
bestimmt. Aufgrund der Form dieser Zerlegungen gilt, dass für jedes
die ersten
Spalten der Matrix
denselben Unterraum aufspannen wie die ersten
Spalten der Matrix
;
zusätzlich sind die Spalten von
orthonormal zueinander. Dieses jedoch ist genau die Situation nach dem
-ten
Schritt einer simultanen Potenziteration. Die Diagonalelemente von
sind wegen
die Approximationen der Eigenwerte von
.
Daher lassen sich die Konvergenzeigenschaften der Potenziteration
übertragen:
Der einfache QR-Algorithmus konvergiert nur, wenn alle Eigenwerte in ihren Beträgen voneinander verschieden sind. Sind die Eigenwerte nach
sortiert, so ist die Konvergenzgeschwindigkeit linear mit einem
Kontraktionsfaktor, der dem Minimum der Quotienten
entspricht.
Insbesondere für reelle Matrizen kann dieser Algorithmus nur konvergieren, wenn alle Eigenwerte reell sind (da sonst konjugiert komplexe Paare mit gleichem Betrag existieren würden). Diese Situation ist für alle symmetrischen Matrizen gegeben.
Der QR-Algorithmus als simultane inverse Potenziteration
Es gilt, falls
invertierbar ist, für die Transponierte
(für komplexe Matrizen die hermitesch Adjungierte) der Inversen
von
und alle ihre Potenzen
Die Inverse einer rechtsoberen Dreiecksmatrix ist wieder eine rechtsobere
Dreiecksmatrix, deren Transponierte eine linksuntere Dreiecksmatrix. Damit
bestimmt der QR-Algorithmus auch eine QL-Zerlegung der Potenzen von .
Aus der Form dieser Zerlegung erhellt, dass für jedes
die letzten
Spalten von
denselben Unterraum aufspannen wie die letzten
Spalten von
.
In der letzten Spalte von
wird somit eine inverse Potenziteration für
durchgeführt, diese Spalte konvergiert also gegen den dualen Eigenvektor zum
kleinsten Eigenwert von
.
Im Produkt
ist also die linke untere Komponente
der sog. Rayleigh-Quotient
der inversen Potenziteration. Die Konvergenzeigenschaften sind analog zum oben
angegebenen.
![Trenner](/button/corpdivider.gif)
![Extern](/button/extern.png)
![Seitenende](/button/stonrul.gif)
© biancahoegel.de
Datum der letzten Änderung: Jena, den: 24.01. 2024