Kapitel 7
Matrizeneigenwertaufgaben
7.1
Aufgabenstellung, Vorbemerkungen
Gegeben sei A ∈ IRN ×N bzw. A ∈ CN ×N . Gesucht sind die Eigenwerte λ ∈ C
und die zugehörigen Eigenvektoren ϕ ∈ CN (ϕ 6= 0), so daß Aϕ = λϕ.
Die Bestimmung von λ und ϕ ist i.a. ein nichtlineares (also infinites) Problem (λ und ϕ werden miteinander multipliziert): G(λ, ϕ) = Aϕ − λϕ = 0. Die
Lösung dieses Problems ist nicht eindeutig bestimmt. Es gibt i.a. mehr als einen
Eigenwert, und jeder Eigenvektor ist nur bis auf eine Konstante bestimmt (man
kann z.B. zusätzlich kϕk = 1 fordern).
Falls ein Eigenwert λ bekannt ist, ist die Bestimmung eines zugehörigen
Eigenvektors ein lineares Problem. Man erhält ihn durch Lösung des unterbestimmten linearen Gleichungssystems (A − λI)ϕ = 0.
Bei bekanntem Eigenvektor ϕ ergibt sich der Eigenwert λ bereits durch
Betrachtung einer einzelnen (skalaren) Gleichung von Aϕ = λϕ.
Bemerkung 7.1 (Charakteristisches Polynom). PA (λ) := det(A − λI) ist
das charakteristische Polynom von A. Es ist ein Polynom N -ten Grades. Die
Nullstellen von PA (λ) sind gerade die Eigenwerte von A: PA (λ) = 0 ⇔ λ ist
Eigenwert von A.
Die Bestimmung der Nullstellen des charakteristischen Polynoms ist eine
Möglichkeit, die Eigenwerte einer Matrix A zu berechnen (z.B. mit dem NewtonVerfahren). Dieser Zugang ist in der Regel jedoch numerisch nicht besonders
interessant.
I
Beim Ansatz von Verfahren zur Bestimmung der Eigenwerte einer Matrix sollte
man die Aufgabenstellung und die eventuell bekannten Eigenschaften der Matrix
berücksichtigen.
Typische Aufgabenstellungen:
1. Gesucht ist ein spezieller Eigenwert (z.B. betraglich größter oder kleinster
Eigenwert) und der zugehörige Eigenvektor (vgl. Abschnitt 7.2).
2. Gesucht sind mehrere Eigenwerte (vgl. Abschnitt 7.3).
3. Gesucht sind alle Eigenwerte (vgl. Abschnitt 7.4, 7.5).
179
180
Matrizeneigenwertaufgaben
In der Praxis typischerweise auftretende Matrizen:
Praktisch auftretende Matrizen sind oft diagonalisierbar, z.B. reell-symmetrisch,
positiv definit o.ä. Hierfür gibt es sehr effiziente Spezialverfahren. Häufig weisen die Matrizen auch eine spezielle Struktur auf und werden zur numerischen
Lösung einer Eigenwertaufgabe bei Differentialgleichungen benötigt. Ein (triviales, numerisch uninteressantes) Beispiel ist:
−u00 = λu
(0 ≤ x ≤ 1),
u(0) = u(1) = 0.
Durch Diskretisierung von −u00 (s. Abschnitt 1.2.4) entsteht eine Matrixeigenwertaufgabe. Typischerweise erhält man dünnbesetzte (“sparse”) Matrizen.
Typische Anwendungsbereiche, in denen Eigenwertaufgaben eine Rolle spielen, sind etwa
– Strukturmechanik (Vibrationen, Resonanz)
– Akustik
– elektrische Netze
– Quantenmechanik (Schrödinger-Gleichung, Eigenzustände)
– Stabilität (dynamischer Systeme, Chaostheorie, Bifurkation) oder
– chemische Reaktionen.
Ein Beispiel einer Brücke (Tacoma Narrows Bridge), die durch Wind in Eigenschwingungen gerät und zusammenstürzt, findet man im Web z.B. unter
http://www.morgenthal.org/index.html.
Bemerkung 7.2 (Diagonalisierbare Matrizen). Wir betrachten spezielle
Typen diagonalisierbarer Matrizen (vgl. Lineare Algebra):
A diagonalisierbar ⇔ Die Jordansche Normalforn besteht nur aus der
Hauptdiagonalen.
⇔ Es gibt eine Basis aus Eigenvektoren (weitere
Hauptvektoren gibt es nicht).
A normal
⇔ A∗ A = AA∗ (A∗ : konjugiert adjungierte Matrix
von A )
A normal
⇒ A diagonalisierbar
A normal
⇔ Es gibt eine Orthonormalbasis aus Eigenvektoren.
Mit anderen Worten: Ist A normal, so ist A “unitär diagonalisierbar”, d.h. es
existiert eine unitäre Matrix U , so daß U ∗ AU eine Diagonalmatrix ist.
A hermitesch (selbstadjungiert) ⇔ A = A∗ (⇒ nur reelle Eigenwerte und
A ist normal)
A reell-symmetrisch
⇔
A = A∗ = AT ∈ IRN ×N
(⇒ nur reelle Eigenwerte)
A unitär
⇔
A∗ = A−1
(⇒ |λ| = 1 für alle Eigenwerte λ). I
Bemerkung 7.3 (Vorbehandlung von Matrizen). Vor der Berechnung der
Eigenwerte einer Matrix überführt man die Matrix typischerweise in eine geeignetere Form. Die Bestimmung der Jordanschen Normalform ist hierbei ungeeignet, da dies ein nichtlineares Problem darstellt. Jede hermitesche Matrix A
ist allerdings (finit, d.h. nach endlich vielen Schritten) durch Ähnlichkeitstransformationen in eine tridiagonale Matrix überführbar, jede quadratische Matrix
181
Matrizeneigenwertaufgaben
ist (finit) in Hessenbergform (vgl. Abb. 7.1) überführbar. Für solche Transformationen sind effiziente Algorithmen verfügbar.
I
0
*
Abbildung 7.1: Hessenbergform einer Matrix
7.2
Vektoriteration und Varianten
Wir betrachten das Problem Aϕ = λϕ. Dabei sei A diagonalisierbar und habe
die Eigenwerte λ1 , . . . , λN mit zugehörigen Eigenvektoren ϕ1 , . . . , ϕN . λ1 sei
“dominierend”, d.h. es gelte
|λ1 | > |λ2 | ≥ |λ3 | ≥ · · · ≥ |λN | ≥ 0.
7.2.1
Vektoriteration
Zu gegebenem x(0) 6= 0 bezeichnet man die Iteration
x(n+1) = Ax(n)
(7.1)
als Vektoriteration (oder Potenzmethode, engl. power method).
Satz 7.1 (Vektoriteration). An x(0) sei ϕ1 beteiligt, d.h. in der Entwicklung
x(0) = c1 ϕ1 + · · · + cN ϕN sei c1 6= 0. Dann gilt:
(i)
x(n) = λ1 n (c1 ϕ1 + r(n) ) mit r(n) → 0 für n → ∞.
³ ´n
³¯ ¯n ´
³ ´n
¯ ¯
λN
λ2
(n)
ϕN = O ¯ λλ21 ¯ .
wobei r = c2 λ1 ϕ2 + · · · + cN λ1
(n+1)
(ii)
x(i)
(n)
x(i)
³¯ ¯n ´
¯ ¯
= λ1 + O ¯ λλ21 ¯
(n+1)
D.h.
x(i)
(n)
x(i)
für alle i = 1, . . . , N mit (ϕ1 )(i) 6= 0.
¯ ¯
¯ ¯
→ λ1 für n → ∞ mit dem Konvergenzfaktor ¯ λλ21 ¯.
Die Vektoriteration liefert also einen zu ϕ1 proportionalen Vektor gemäß (i)
und den Eigenwert λ1 gemäß (ii).
Beweis: zu (i):
x(n) = An x(0) = An (c1 ϕ1 + . . . + cN ϕN )
= λn1 c1 ϕ1 + . . . + λnN cN ϕN = λn1 (c1 ϕ + r(n) ).
182
Matrizeneigenwertaufgaben
zu (ii):
(n+1)
x(i)
(n)
x(i)
³
´P
³ ´n
³³
´
´
PN ³ λj ´n
λj
λj
N
λn+1
c
ϕ
+
c
ϕ
−
1
c
ϕ
+
1
1
j
j
j
j
1
j=2 λ1
j=2 λ1
λ1
(i)
(i)
´
³
³ ´n
=
P
λj
N
n
λ1 c1 ϕ1 + j=2 λ1 cj ϕj
(i)
µ
µ¯ ¯n ¶¶
¯ λ2 ¯
.
= λ1 1 + O ¯¯ ¯¯
λ1
¥
Beachte: Die Forderung c1 6= 0 ist numerisch unerheblich, denn falls c1 = 0 sein
sollte, bewirken Rundungsfehler i.d.R., dass c1 6= 0 bei x(1) , x(2) , etc.
7.2.2
Skalierungen
Die dauernde Veränderung der Größenordnung von x(n) um den asymptotischen
Faktor λ1 kann durch geeignete Skalierungen vermieden werden.
Durch y (n) := x(n)1 x(n) erhält man eine Skalierung der i-ten Komponente
( )(i)
(n)
von x
auf 1 (im Falle (ϕ1 )(i) 6= 0 ist für alle hinreichend großen n auch
©
¡ (n) ¢
ª
6= 0). Wir bilden daher die Folge y (n) gemäß
x
(i)
y (n+1) =
=
mit µ(n+1) =
¡
¡
1
Ax(n)
1
Ay (n)
1
(Ay(n) )(i)
¢
(i)
¢
Ax(n) = ¡
Ay (n) =
(i)
und y (0) =
x(n)
1
³
´
1
(n)
¡
¢
x
Ax(n)
(n)
(i)
Ax
(i)
(i)
¢
µ(n+1)
Ay (n)
x(0)
.
x
( (0) )(i)
Hierfür lauten die Aussagen des Satzes 7.1
y
(n)
µ(n)
µ¯ ¯n ¶
¯ λ2 ¯
1
=
ϕ1 + O ¯¯ ¯¯
(ϕ1 )(i)
λ1
ï ¯
!
¯ λ2 ¯n−1
= λ1 + O ¯¯ ¯¯
.
λ1
In der Praxis am günstigsten ist die Division durch die betraglich größte Komponente.
Allgemeine Skalierungsmöglichkeiten:
(n)
(n+1)
,`
i
Andere Skalierungen erhält man, wenn man in (7.2) µ(n+1) = hAx
mit
hx(n) ,`(n) i
©
ª
einem konstanten Vektor `(n+1) = ` oder einer Folge `(n) wählt. Der Wahl
in (7.2) entspricht ` = ei , denn
¡ (n) ¢
³
´
Ax
hAx(n) , ei i
(i)
(n)
=
Ay
= ¡ (n) ¢
.
(i)
hx(n) , ei i
x
(i)
Ein weiterer gebräuchlicher Spezialfall ist `(n+1) = ` = x(0) :
183
Matrizeneigenwertaufgaben
A sei hermitesch und ϕ1 , . . . , ϕN sei eine Orthonormalbasis aus Eigenvektoren.
Dann ergibt sich
µ(n+1) =
a(n+1)
hAx(n) , x(0) i
=
= q (n)
hx(n) , x(0) i
a(n)
mit den Schwarzschen Konstanten a(n) und den Schwarzschen Quotienten q (n) .
Da A hermitesch ist, gilt
a(n) = hAx(n−1) , x(0) i = hAk x(0) , Am x(0) i mit
k + m = n.
Man beachte hierbei, dass zur Bestimmung von q (2n) nur x(n) benötigt wird:
q (2n) =
a(2n+1)
hAx(n) , x(n) i
=
=: R(n)
a(2n)
hx(n) , x(n) i
mit dem Rayleighquotienten R(n) . Analog zu Satz 7.1 (ii) läßt sich zeigen
ï ¯ !
¯ λ2 ¯2n
(n)
.
R = λ1 + O ¯¯ ¯¯
λ1
¯ ¯2
¯ ¯
Man hat also jetzt einen Konvergenzfaktor von ¯ λλ12 ¯ .
Als dritten Spezialfall betrachten wir noch `(n+1) = x(n+1) :
A sei wieder hermitesch und ϕ1 , . . . , ϕN eine Orthonormalbasis aus Eigenvektoren. Dann erhält man
hAx(n) , Ax(n) i
a(2n+2)
=
= q (2n+2) q (2n+1)
hx(n) , x(n) i
a(2n)
³¯ ¯´
¯ ¯
= λ21 + O ¯ λλ21 ¯ .
µ(n+1) =
und damit µ(n)
7.2.3
Shiftstrategien
¯ ¯
¯ ¯
Wir wollen uns nun mit der Frage beschäftigen, wie der Konvergenzfaktor ¯ λλ12 ¯
verbessert werden kann. Hierzu ist folgende einfache Überlegung hilfreich:
Aϕ = λϕ
⇔ (A − αI)ϕ = (λ − α)ϕ
↑
↑
Eigenwerte von A
Eigenwerte von A − αI
λ1 , . . . , λN
λ1 − α, . . . , λN − α
Konvergenzbeschleunigung erreicht man hiermit leicht, falls z.B. A hermitesch und positiv definit ist und falls gilt
λ1 > λ2 ≥ λ3 ≥ . . . ≥ λN > 0.
Dann bestimmt man α so, daß
¯
¯ λ2 − α
¯
¯ λ1 − α
¯
¯
¯
¯
¯ < ¯ λ2
¯ λ1
¯
¯
¯
¯.
¯
184
Matrizeneigenwertaufgaben
(n+1) = (A − αI)x(n) konvergiert dann mit dem KonvergenzDas Verfahren
¯ x
¯ ¯
¯
¯ −α ¯
(n+1) = Ax(n) (Konvergenzfaktor ¯ λ2 ¯).
und
damit
schneller
als
x
faktor ¯ λλ12 −α
¯
¯ λ1 ¯
Konvergenzerzeugung ist z.B. im Fall
|λ1 | = λ1 = −λ2 > |λ3 | ≥ . . .
möglich, in dem keine Dominanz von λ1 vorliegt. Durch eine geeignete Wahl
von α (α ≈ λ2 ) wird λ1 − α dominierend.
7.2.4
Inverse Iteration
Für die Berechnung des betraglich kleinsten Eigenwertes setzen wir voraus, dass
|λ1 | ≥ |λ2 | ≥ . . . ≥ |λN −1 | > |λN | > 0. Dann gelten für die inverse Iteration
Ax(n+1) = x(n)
(7.2)
alle obigen Aussagen und Überlegungen für λN statt λ1 , denn λN ist der betragsgrößte Eigenwert von A−1 . Praktisch wird man A nicht invertieren, sondern
in jedem Schritt das Gleichungssystem (7.2) lösen (etwa mit fester LR-Zerlegung
und sich ändernden rechten Seiten x(n) ). Insbesondere kann A−1 für dünnbesetzte Matrizen A vollbesetzt sein, so daß der Aufwand über die Berechnung
der Inversen viel zu hoch wäre.
Die Idee der Verschiebung ist auch bei der inversen Vektoriteration anwendek bekannt (z.B. berechnet). Dann
bar (Wieland). Für λk sei eine Näherung λ
e
wählt man α = λk und iteriert gemäß
(A − αI)x(n+1) = x(n) .
ek ≈ λk ist 1 (sehr)
Ein Eigenwert von (A−αI)−1 ist dann λk1−α . Wegen α = λ
λk −α
groß und daher in der Regel stark dominierend. Man darf also sehr günstige
Konvergenz erwarten und kann so schnell eine verbesserte Näherung für λk
berechnen.
7.3
Deflation, Treppeniteration
In diesem Abschnitt betrachten wir den Fall, dass wir nicht nur einen sondern
mehrere Eigenwerte berechnen wollen.
Wir wollen zunächst λ2 berechnen. Hierzu sei A hermitesch und ferner gelte
|λ1 | > |λ2 | > |λ3 | ≥ · · · ≥ |λN | ≥ 0.
Berechne zunächst λ1 , ϕ1 und bilde dann
x
e(n+1) := Ax(n)
x(n+1) = x
e(n+1) − he
x(n+1) , ϕ1 iϕ1 ,
(7.3)
das heißt x(n+1) ist orthogonal zu ϕ1 , also ist ϕ1 nicht an x(n+1) beteiligt.
185
Matrizeneigenwertaufgaben
Wie sieht jetzt die Konvergenz gegen ϕ2 bzw. λ2 aus? Wir können schreiben
x(n+1) = Ax(n) − hAx(n) , ϕ1 iϕ1 = Ax(n) − λ1 hx(n) , ϕ1 iϕ1
= (A − λ1 ϕ1 ϕ1 T ) x(n) .
|
{z
}
=:W
Die Matrix W ist ebenfalls hermitesch. Offenbar gilt
W ϕ1 = 0,
W |[ϕ2 ,...,ϕN ] = A|[ϕ2 ,...,ϕN ] ,
das heißt W hat die Eigenwerte 0, λ2 , . . . , λN . Die Vektoritertion mit der
¯
¯ Matrix
¯ ¯
W liefert also λ2 und ϕ2 gemäß Satz 7.1 mit dem Konvergenzfaktor ¯ λλ32 ¯.
Beim obigen Vorgehen wird erst ϕ1 , λ1 berechnet und dann ϕ2 , λ2 . Besser ist
es, beide gleichzeitig zu berechnen. Dazu setzen wir in (7.3) statt ϕ1 selbst die
jeweils bekannte Näherung für ϕ1 ein. Mit Skalierung ergibt sich dann folgendes
Verfahren:
y (n+1)
=
w(n+1) =
1
µ(n+1)
1
σ (n+1)
Ay (n)
¡
(→ ϕ1 , λ1 )
Aw(n) − hAw(n) , y (n+1) iy (n+1)
¢
(→ ϕ2 , λ2 ).
Bemerkung 7.4 (Treppeniteration). Die Fortführung obiger Idee wird als
Treppeniteration bezeichnet. Falls A hermitesch und |λ1 | > |λ2 | > · · · > |λN | >
0, dann lassen sich so alle Eigenwerte berechnen.
I
7.4
LR- und QR-Algorithmus
LR- und QR-Algorithmen sind die am meisten verwendeten Algorithmen zur
Berechnung sämtlicher Eigenwerte einer Matrix A. Jedoch werden sie in der
Regel nicht direkt auf beliebige Matrizen A angewendet, sondern die Matrix
wird zunächst vorbehandelt und in eine Form gebracht, die für Konvergenz und
Algorithmik vorteilhaft ist.
7.4.1
Der LR-Algorithmus
Der LR-Algorithmus wird hier beschrieben, ohne genaue Voraussetzungen anzugeben. Implizit wird angenommen, daß der Algorithmus (auch ohne Pivotwahl)
durchführbar ist.
Ist
A = A0 = L0 R0
die LR-Zerlegung der Matrix A, so bilden wir
A1 := R0 L0
durch einfache Matrixmultiplikation und führen für A1 wieder eine LR-Zerlegung
A1 = L1 R1
186
Matrizeneigenwertaufgaben
durch usw. Allgemein erhalten wir damit
An = Ln Rn
An+1 = Rn Ln
Satz 7.2 (LR-Algorithmus). Es sei A ∈ CN ×N . Der LR-Algorithmus sei
durchführbar. Für die Eigenwerte von A gelte
|λ1 | > |λ2 | > . . . > |λN |.
(7.4)
Für die Matrizen φ, φ−1 , durch die A auf die Jordansche Normalform
λ1
0
−1
..
D=
= φ Aφ
.
0
λN
transformiert wird, gebe es Dreieckszerlegungen
φ = Lφ Rφ ,
φ−1 = Lφ−1 Rφ−1 .
Dann konvergieren die Matrizen An , Ln und Rn :
λ1
∗
..
lim An = lim Rn =
,
.
n→∞
n→∞
0
λN
lim Ln = I.
n→∞
¯ ¯
¯ ¯
Die Konvergenz ist linear und hängt vom Verhalten der Quotienten ¯ λλki ¯ ab.
Beweis: vgl. Stoer-Bulirsch.
¥
Bemerkung 7.5 (Bezug zur Treppeniteration). Der LR-Algorithmus ist
theoretisch äquivalent zur Treppeniteration.
I
Bemerkung 7.6 (Erläuterungen zu Satz 7.2).
(1) Die Voraussetzung (7.4) ist abschwächbar.
(2) Falls keine LR-Zerlegung möglich ist, dann kann man u.U. mit Pivotisierung weiterkommen.
(3) Bei symmetrisch positiv definiten Matrizen ist statt der LR-Zerlegung
auch eine entsprechende Cholesky-Variante möglich und sinnvoll (weniger
Rechenaufwand).
(4) Eine Konvergenzbeschleunigung ist durch Verschiebungen A → A +
− αI
( Shift-Techniken“) erreichbar.
”
(5) Der LR-Algorithmus ist für vollbesetzte Matrizen sehr rechenaufwändig.
Daher ist es ratsam, den LR-Algorithmus nicht auf die ursprüngliche Matrix A anzuwenden, sondern diese zunächst vorzubehandeln (z.B. diese
mit den Verfahren von Givens (vgl. Abschnitt 7.5.2) oder Householder
auf Tridiagonal- oder Hessenbergform zu transformieren).
(6) Die numerische Problematik (Stabilität) von LR-Zerlegungen schlägt sich
auch beim LR-Algorithmus nieder.
I
187
Matrizeneigenwertaufgaben
Der QR-Algorithmus
7.4.2
Eine Weiterentwicklung des LR-Algorithmus ist der QR-Algorithmus, der auf
orthogonalen Zerlegungen einer Matrix beruht.
Der QR-Algorithmus zur Eigenwertberechnung lautet
A0 = A
An = Qn Rn
An+1 = Rn Qn
(7.5)
(n = 0, 1, 2, . . .)
mit
Q∗n Qn = I
∗
..
.
∗
∗ ···
..
R=
.
(d.h. Qn ist unitär) und
0
Beim QR-Verfahren verwenden wie also eine QR-Zerlegung statt einer LRZerlegung. Die Berechnung der QR-Zerlegung erfolgt mit Householder-Transformationen wie in Abschnitt 6.2 erläutert.
Satz 7.3 Die im QR-Verfahren auftretenden Matrizen An sind alle ähnlich zur
Matrix A.
Beweis:
Qn An+1 Q∗n = Qn Rn Qn Q∗n = Qn Rn = An .
(7.6)
Eine wiederholte Anwendung obiger Beziehung liefert die Behauptung.
¥.
Den folgenden Satz beweisen wir für den Fall reell symmetrischer Matrizen.
Ähnliche, verallgemeinerte Aussagen sind auch für komplexe unsymmetrische
Matrizen möglich.
Satz 7.4 (QR-Algorithmus). Sei A ∈ IRN ×N symmetrisch. A habe N Eigenwerte λ1 , . . . , λN mit
|λ1 | > |λ2 | > · · · > |λN | > 0.
(7.7)
³
(n)
Die An , Qn , und Rn seien wie in (7.5). Dann gilt mit An = aij
(i)
lim Qn = I,
n→∞
(ii)
lim Rn = Λ = diag (λ1 , . . . , λN ),
n→∞
(iii)
(n)
aij
µ¯ ¯n ¶
¯ λi ¯
= O ¯¯ ¯¯
λj
für
i > j.
´
i,j=1,...,N
:
188
Matrizeneigenwertaufgaben
Beweis: Wir zeigen zunächst per vollständiger Induktion die Beziehung
An = Q1 · · · Qn Rn · · · R1
| {z } | {z }
n = 1, 2, . . .
(7.8)
=:Un
=:Pn
Der Fall n = 1 ist klar wegen A = A1 = Q1 R1 .
Außerdem gilt
An+1 =
=
=
=
=
Qn+1 Rn+1
QTn An Qn
QTn QTn−1 An−1 Qn−1 Qn
QTn · · · QT1 AQ1 · · · Qn
Pn −1 APn
(wegen (7.6))
(wieder wegen (7.6))
(wieder wegen (7.6))
(wegen (7.8))
und deshalb APn = Pn An+1 .
Mit dieser Beziehung können wir jetzt den Induktionsschritt durchführen.
Es ist
An+1 = AAn = APn Un = Pn An+1 Un = Pn Qn+1 Rn+1 Un = Pn+1 Un+1 .
(7.9)
Damit haben wir die Beziehung (7.8) bewiesen. Weil Pn orthogonal und Un
rechte obere Dreiecksmatrix ist, bedeutet dies, dass wir die QR-Zerlegung von
An = Pn Un von An durch die QR-Zerlegungen der A1 , . . . , An ausdrücken
können. Weiter gilt
An = QΛn QT
mit Λn = diag(λn1 , . . . , λnN ).
Wir nehmen nun der Einfachheit halber an, QT hätte eine LR-Zerlegung mit
diag(L) = (1) (sonst können wir dies erreichen, indem wir A geeignet permutieren). Dann ist
An = QΛn LR = Q (Λn LΛ−n )(Λn R)
| {z }
³ ´n
(Λn LΛ−n )ij = lij λλji
→ 0 für
n → ∞, j > i,
(7.10)
denn L ist untere Dreiecksmatrix, und für i > j (die Nichtdiagonalelemente)
ist nach Voraussetzung λj > λi . Damit haben wir
Λn LΛ−n = I + En
mit En → 0 für n → ∞.
Einsetzen in (7.10) liefert
An = Q(I + En )Λn R.
Nun wenden wir auch auf In + En (formal) eine QR-Zerlegung
en R
en
I + En = Q
en positiv seien. Aus der Eindeutigkeit der
an, wobei die Diagonalelemente von R
QR-Zerlegung und limn→∞ En = 0 folgt
en , R
en → I
Q
für n → ∞.
189
Matrizeneigenwertaufgaben
Damit haben wir eine zweite QR-Zerlegung von An hergeleitet, da
e n )(R
en Λn R).
An = (QQ
Bis auf die Vorzeichen in der Diagonalen gilt daher (Vergleich mit (7.8))
e n , Un = R
en Λn R.
Pn = QQ
Für n → ∞ folgt
T
Qn = QTn−1 · · · QT1 Q1 · · · Qn−1 Qn = Pn−1
Pn
| {z }
=I
|
{z
}
=I
e Tn−1 QT QQ
en = Q
e Tn−1 Q
en → I
= Q
−1
⇒ (i),
−1
Rn = Rn · · · R1 R1 · · · Rn−1 = Un Un−1 −1
en Λn RR−1 Λ−(n−1) R
e−1 = R
e n ΛR
e−1 → Λ
= R
n−1
n−1
⇒ (ii).
Schließlich ergibt sich (iii) aus
lim An = lim Qn Rn = lim Rn = Λ
n→∞
n→∞
n→∞
zusammen mit (7.10).
¥
Bemerkung 7.7 Es gilt ferner unter den gleichen Voraussetzungen wie in
Satz 7.4:
(1) A symmetrisch ⇒ An symmetrisch.
(2) A symmetrisch und tridiagonal ⇒ An symmetrisch und tridiagonal.
Einen Beweis hierfür findet man z.B. in Deuflhard/Hohmann.
I
Bemerkung 7.8 (Erläuterungen zu Satz 7.4).
(1) Die Voraussetzung (7.4) ist abschwächbar.
(2) Eine geeignete Shiftstrategie (mit mehreren Shiftparametern) liefert quadratische Konvergenz.
(3) Der Rechenaufwand ist bei der QR-Zerlegung pro Schritt etwa doppelt so
groß wie bei der LR-Zerlegung.
(4) Die Konvergenz des QR-Verfahrens ist im Normalfall doppelt so schnell
wie die des LR-Verfahrens. Die numerischen Eigenschaften des QR-Verfahrens sind in der Regel besser.
(5) Auch die QR-Zerlegung sollte man nicht naiv anwenden, sondern erst
nach Vorbehandlung von A, indem man A zunächst in Tridiagonal- oder
Hessenberggestalt bringt. Der Aufwand ist sonst zu hoch.
I
190
Matrizeneigenwertaufgaben
7.5
Die Verfahren von Jacobi und Givens
A ∈ IRN ×N sei symmetrisch. Sowohl das Verfahren von Jacobi als auch das
Verfahren von Givens beruht auf orthogonalen Transformationen mithilfe von
Givensrotationen, vgl. Abschnitt 6.2.3.
7.5.1
Die Idee des Verfahrens von Jacobi
Das Verfahren von Jacobi besteht darin, Givensrotationen durchzuführen, z.B.
in der Reihenfolge
Ω12 , Ω13 , · · ·
Ω1N
Ω23 , · · ·
Ω2N
..
..
.
.
ΩN −1,N
mit zyklischer Wiederholung. Letztere ist notwendig, weil die in den vorangegangenen Schritten erzeugten Nullen nicht erhalten bleiben. Bei zyklischer
Wiederholung der Schritte gilt
λ1
0
..
A(n) → D =
.
.
0
λN
Die Quadratsumme der Nichtdiagonalelemente
X (n)
S(A(n) ) =
|ak` |2
k6=`
ist ein Maß für die Konvergenz dieses Verfahrens. Aus den Formeln für die
Einträge c, s der Givensrotationen (s. (6.3)) folgert man
(n)
0 ≤ S(A(n+1) ) = S(A(n) ) − 2|ak` |2 ≤ S(A(n) ).
Wir haben also eine monoton fallende Folge. Ohne Beweis vermerken wir:
S(A(n) ) → 0.
Unter besonderen Umständen (A hat nur einfache Eigenwerte) erhält man sogar quadratische Konvergenz. Eine Variante dieses Verfahrens besteht darin, die
Schritte nicht in einer festen Reihenfolge mit zyklischer Wiederholung durchzuführen, sondern stattdessen jeweils das betragsgrößte Element ak` mit k 6= `
zu Null zu machen. Dieser Ansatz führt zu einer Konvergenzbeschleunigung.
7.5.2
Die Idee des Verfahrens von Givens
Das Verfahren von Givens ist ein finites Verfahren. Es beruht auf der gleichen
Idee wie das Jacobi-Verfahren, aber hier sorgt man dafür, daß einmal erzeugte
Nullen nicht wieder ungleich Null werden. Damit besteht das Verfahren nur aus
endlich vielen Schritten. Bei diesem Verfahren gelangt man allerdings nicht zu
191
Matrizeneigenwertaufgaben
einer Diagonalmatrix, sondern nur zu einer Tridiagonalmatrix. Die Reihenfolge
der Rotationen bei diesem Verfahren ist
Ω23 , Ω24 , · · ·
Ω34 , · · ·
..
.
Ω2N
Ω3N
..
.
.
ΩN −1,N
Allerdings wird hier durch die Transformation mit Ωk` nicht das Element ak` zu
Null transformiert, sondern (durch geeignete Wahl von θ) das Element a`,k−1
und damit auch ak−1,` (im reell-symmetrischen Fall).
Eine Verfeinerung des Verfahrens von Givens ist das Householder-Verfahren
(vgl. Stoer-Bulirsch).
Beispiel 7.1 Als Beispiel betrachten wir
a11 a12
a21 a22
A=
a31 a32
a41 a42
Die Reihenfolge der Rotationen ist dann
den Fall einer 4 × 4-Matrix
a13 a14
a23 a24
.
a33 a34
a43 a44
Ω23 , Ω24
.
Ω34
Durch die Transformation mit Ω23 werden die Matrixelemente a31 und a13 auf
Null gebracht, Ω24 führt dann auf eine Matrix der Gestalt
∗ ∗ 0 0
∗ ∗ ∗ ∗
0 ∗ ∗ ∗ .
0 ∗ ∗ ∗
Ω34 bringt dies dann auf Tridiagonalgestalt
∗ ∗ 0 0
∗ ∗ ∗ 0
0 ∗ ∗ ∗
0 0 ∗ ∗
.
Allgemein entsteht bei diesem Verfahren im Sinne der obigen Reihenfolge eine
−2)
Schritten.
J
Tridiagonalmatrix nach (N −1)(N
2
7.6
Eigenwertabschätzungen
Eigenwertabschätzungen können theoretisch sehr ausführlich behandelt werden
(vgl. z.B. Stoer/Bulirsch). Sie sind in vielen Feldern, darunter Funktionalanalysis und Quantenmechanik von grundsätzlicher Bedeutung. Auch wir haben im
192
Matrizeneigenwertaufgaben
Zusammenhang mit dem Spektralradius bereits einige Eigenwertabschätzungen
kennengelernt. So war etwa für eine Iterationsmatrix M
p
ρ(M ) = lim n kM n k
n→∞
= inf {kM k : k · k passend}
ρ(M ) ≤ kM k ,
7.6.1
oder
k · k passend.
Der klassische Satz von Gerschgorin
Die Diagonalelemente einer Matrix sind Näherungen für die Eigenwerte, und
zwar umso besser, je stärker diagonaldominant die Matrix ist. Dies ist die zentrale Aussage des Satzes von Gerschgorin. Genau lautet er:
Satz 7.5 (Satz von Gerschgorin). Es sei A ∈ CN ×N . Bilden wir die Kreise
von Gerschgorin
Zi := {µ ∈ C : |µ − aii | ≤ zi }
Sj := {µ ∈ C : |µ − ajj | ≤ sj }
mit
zi
=
X0
k
sj
=
X0
N
X
|aik |,
|aik | =
k=1
k6=i
|akj | =
k
N
X
|akj |,
k=1
k6=j
dann gilt:
(1) Jeder Eigenwert von A liegt in Z :=
SN
i=1 Zi .
(2) Jede Zusammenhangskomponente von Z enthält so viele Nullstellen des
charakteristischen Polynoms von A, wie Kreise Zi an ihr beteiligt sind.
S
(3) Diese Aussagen gelten entsprechend für Sj , S := N
j=1 Sj .
Beweis: Wir zeigen hier nur (1); für (2) und (3) vgl. etwa Stoer/Bulirsch.
Sei λ Eigenwert von A und k · k = k · k∞ . Sei ϕ ein zugehöriger Eigenvektor
zum Eigenwert λ. D bezeichne die “Diagonalmatrix” von A. Ferner sei λ 6= aii
(für alle i). Dann ist
(A − D)ϕ = (λI − D)ϕ mit (λI − D) regulär wegen λ 6= aii .
Es gilt
(λI − D)−1 (A − D)ϕ = ϕ.
Der Übergang zur Norm liefert
k(λI − D)−1 (A − D)k∞ kϕk∞ ≥ kϕk∞
(mit k · k passend), d.h.
k(λI − D)−1 (A − D)k∞ ≥ 1.
193
Matrizeneigenwertaufgaben
Hieraus folgt sofort
max
1≤i≤N
P0
k |aik |
≥ 1.
|λ − aii |
Da λ Eigenwert ist, existiert ein i, so daß
|λ − aii | ≤
0
X
|aik |.
k
Dies bedeutet, λ liegt in einem der Gerschgorin-Kreise, und hieraus folgt die
Aussage (1) des Satzes.
Die Aussage (2) folgt mit Stetigkeitsargumenten. Man betrachtet dazu die
Schar von Matrizen
At = D + t(A − D) mit
At=0 = D
und At=1 = A
und benutzt die Tatsache, daß die Eigenwerte stetig von t abhängen (vgl. Stoer/Bulirsch II).
¥
Bemerkung 7.9 Im Falle lauter disjunkter Kreise liefert der Satz brauchbare
Abschätzungen. Diese sind z.B. anwendbar, wenn A bereits eine gute Näherung
für eine Diagonalmatrix ist. Man kann diese Aussagen auch als Abbruchkriterium verwenden (z.B. im Jacobi-Verfahren).
I
7.6.2
Rayleigh-Quotienten
Satz 7.6 (ohne Beweis).
A ∈ CN ×N sei hermitesch. Die (reellen) Eigenwerte von A seien so angeordnet,
daß
λ1 ≥ λ 2 ≥ · · · ≥ λ N .
Dann ist
hAx, xi
x6=0 hx, xi
hAx, xi
.
=
min
x∈CN , x6=0 hx, xi
λ1 =
λN
max
x∈CN ,
und
Allgemeiner gilt sogar
λk+1 =
7.7
min
y1 ,...,yk
∈CN
max
x∈CN , x6=0
(x,yκ )=0 für κ≤k
hAx, xi
.
hx, xi
Rekapitulation
Bzgl. eines zusammenfassenden Überblicks auf Verfahren zu Eigenwertberechnungen verweisen wir auf Abbildung 7.2. Insgesamt können wir die folgenden
prinzipiellen Zugänge zu Eigenwertberechnungen unterscheiden:
194
Matrizeneigenwertaufgaben
A:
H:
T:
D:
D’:
D'
Ausgangsmatrix (allgemeine quadratische Matrix)
Hessenbergform
Tridiagonalform
Dreiecksmatrix
Diagonalmatrix
A
< He >
(LR)
(QR)
H
LR
QR
La
LR
QR
D
La
: Verf. setzt A als hermitesch voraus
> : finites Verfahren
Ho :
He :
La :
Gi :
Ja :
Verfahren
Verfahren
Verfahren
Verfahren
Verfahren
von
von
von
von
von
<
T
Householder (vgl. Stoer-Bulirsch)
Hessenberg
Lanczos
Givens
Jacobi
Abbildung 7.2: Überblick über klassische Verfahren zur Bestimmung von Eigenwerten einer Matrix A.
• Berechung über das charakteristische Polynom (i.a. völlig ineffizient, kann
aber für vorbehandelte Matrizen, die in Tridiagonal- oder Hessenbergform
gebracht wurden, durchaus von Interesse sein),
• Berechnung über Ähnlichkeitstransformationen,
• Berechnung über Lösung des nichtlinearen Gleichungssystems
G(λ, ϕ) = Aϕ − λϕ = 0 ,
kϕk = 1
(das Jacobi-Verfahren läßt sich z.B. auch als ein solches Lösungsverfahren
deuten),
• Berechnung über die Vektoriteration (Potenzmethode) und Varianten.
© Copyright 2024 ExpyDoc
