Numerische Lineare Algebra - tu-chemnitz.de · Jacobi-Verfahren Konvergenz,zyklischeIndexwahl...

Numerische Lineare Algebra

Oliver Ernst

Professur Numerische Mathematik

Wintersemester 2013/14

Inhalt I

1 Organisatorisches2 Einleitung2.1 Vorbemerkungen2.2 Eigenwertaufgaben in den Anwendungen3 Das QR-Verfahren3.1 Vektoriteration3.2 Unterraumiteration3.3 Reduktion auf Hessenberg-Gestalt3.4 QR Iteration mit impliziten Shifts3.5 Konvergenz der Unterraumiteration3.6 Konvergenz des QR-Verfahrens4 Hermitesche Eigenwertaufgaben4.1 Jacobi-Verfahren4.2 Bisektion4.3 Divide & Conquer - Verfahren4.4 Berechnung der Singulärwertzerlegung5 Anhang: Grundlagen aus der Linearen Algebra5.1 EigenwerteOliver Ernst (Numerische Mathematik) Numerische Lineare Algebra Wintersemester 2013/14 2 / 177

Inhalt II

5.2 Störungstheorie5.3 Householder-Transformationen und Givens-Rotationen5.4 Winkel und Abstand zwischen Unterräumen

Oliver Ernst (Numerische Mathematik) Numerische Lineare Algebra Wintersemester 2013/14 3 / 177

Inhalt

1 Organisatorisches

2 Einleitung

3 Das QR-Verfahren

4 Hermitesche Eigenwertaufgaben

5 Anhang: Grundlagen aus der Linearen Algebra


Hermitesche EigenwertaufgabenEinleitung

Für Hermitesche Matrizen A = AH ∈ Cn×n ergeben sich, wie wir in unserenbisherigen Untersuchungen schon festgestellt haben, einige Vereinfachungen bei derBerechnung von Eigenwerten:

Reelle Eigenwerte, orthogonale Eigenvektoren, diagonalisierbar,Eigenwerte/-vektoren stabiler gegen Störungen von A,Interlacing-Eigenschaft der Eigenwerte,QR: aus Hessenberg wird Hermitesch tridiagonal,QR: aus quadratischer Konvergenz wird kubische Konvergenz.

Neben dem QR-Verfahren existieren im Hermiteschen Fall einige weitere Verfah-ren, die in manchen Situationen sogar dem QR-Verfahren vorzuziehen sind. Wirbetrachten hier den reellen (symmetrischen) Fall.


Inhalt

1 Organisatorisches

2 Einleitung

3 Das QR-Verfahren

4 Hermitesche Eigenwertaufgaben4.1 Jacobi-Verfahren4.2 Bisektion4.3 Divide & Conquer - Verfahren4.4 Berechnung der Singulärwertzerlegung



Jacobi-VerfahrenGrundidee

Bei diesem auf C.G.J. Jacobi im Jahr 1845 zurückgehenden Verfahren wird A uni-tären Ähnlichkeitstransformationen unterzogen mit dem Ziel, alle Außendiagonal-einträge hinreichend zu verkleinern. Als Maß für deren Größe führen wir ein

off(A) := ‖A− diag(A)‖F =

n∑j=1

n∑k=1k 6=j

|aj,k|2

1/2

Die orthogonalen Transformationen beruhen (vgl. Givens-Rotationen) auf Rotatio-nen in 2 achsenparallelen Unterräumen: im 2× 2-Fall lässt sich[

c s−s c

]T [a bb d

] [c s−s c

]=

[λ1 00 λ2

], 0 ≤ s, c ≤ 1, s2 + c2 = 1,

erreichen durch Wahl von c = 1/√

1 + t2, s = tc mit t die betragskleinere Nullstellevon

t2 + 2τt− 1 = 0, t =s

c, τ :=

aqq − app2apq

.


Jacobi-VerfahrenRotationsmatrix

Das Durchführen der (c, s)-Rotation, c = cosϑ, s = sinϑ, in der (p, q)-Koordinaten-ebene wird bewerkstelligt durch

J = J (p, q, ϑ) :=

1 · · · 0 · · · 0 · · · 0...

. . ....

......

0 · · · c · · · s · · · 0...

.... . .

......

0 · · · −s · · · c · · · 0...

......

. . ....

0 · · · 0 · · · 0 · · · 1

p

q

p q

und damit den Übergang A 7→ B := JTAJ mit J = J (p, q, ϑ). Die neue MatrixB unterscheidet sich von A nur in den beiden Zeilen und Spalten mit Indices p undq.


Jacobi-VerfahrenReduktion der Nebendiagonalen

Da orthogonale Transformationen die Frobenius-Norm invariant lassen, gilt wegen[bpp bpqbqp bqq

]=

[c s−s c

]T [app apqaqp aqq

] [c s−s c

]zunächst

a2pp + a2qq + 2apq = b2pp + b2qq + 2bpq = b2pp + b2qq

und damit auch

off(B)2 = ‖B‖2F −n∑

i=1

b2ii = ‖A‖2F −n∑

i=1

b2ii

= ‖A‖2F −n∑

i=1

a2ii + (a2pp + a2qq − b2pp − b2qq)

= off(A)2 − 2a2pq.

In diesem Sinne rückt A mit jedem Jacobi-Schritt näher an Diagonalform.


Jacobi-VerfahrenIndexwahl

Es ist

‖B −A‖2F = 4(1− c)n∑

i=1i6=p,q

(a2ip + a2iq) +2

c2a2pq.

Die Wahl der kleineren Nullstelle stellt |ϑ| ≤ π/4 sicher und minimiert den obigenAbstand zwischen A und B .

Sind die Rotationsparanmeter p, q und ϑ bestimmt, erfordert die Aufdatierung

A← J (p, q, ϑ)TAJ (p, q, ϑ)

lediglich 6n arithmetische Operationen.

Im klassischen Jacobi-Verfahren wählt man in jedem Schritt die Indices p, q, umoff(A) möglichst stark zu reduzieren, d.h. so, dass |apq| maximal mit p 6= q.

Der folgende Algorithmus überschreibt A = AT ∈ Rn×n mit VTAV , wobei V dieakkumulierte orthogonale Transformation enthält.


Jacobi-VerfahrenAlgorithmus, lineare Konvergenz

Algorithmus 3 : Klassisches Jacobi-Verfahren.Gegeben : A = AT ∈ Rn×n, tol > 0.

1 V ← In; ε← tol ‖A‖F2 while off(A) > ε do3 Wähle p, q so, dass |apq| = maxi 6=j |aij |4 Bestimme (c, s) um apq im 2× 2-Fall auszunullen5 A← J (p, q, ϑ)TAJ (p, q, ϑ)6 V ← VJ (p, q, ϑ)

Satz 4.1Bezeichnet Ak die k-te Iterierte des klassischen Jacobi-Verfahrens, so gilt

off(Ak)2 ≤(

1− 1

N

)k

off(A)2, N =n(n− 1)

2.

In diesem Sinne liegt lineare Konvergenz vor.


Jacobi-VerfahrenKonvergenz, zyklische Indexwahl

Weitergehende Resultate von [Schönhage, 1964] sowie [van Kempen, 1966] zeigen,dass das klassische Jacobi-Verfahren sogar quadratisch konvergiert.

Die Anzahl der zu einer vorgegebenen Reduktion von off(A) erforderlichen Schrittedes klassischen Jacobi-Verfahrens geben [Brent & Luk, 1985] heuristisch als propor-tional zu log n an.

Zyklisches Jacobi-Verfahren: Das Missverhältnis zwischen O(n) Rechenoperationenpro Jacobi-Schritt gegenüber O(n2) Aufwand bei der Suche nach dem maximalen|apq| wird beim zyklischen Jacobi-Verfahren dadurch vermieden, dass man zyklischalle Nebendiagonalelemente durchläuft, z.B. in der Reihenfolge (n = 4)

(1, 2), (1, 3), (1, 4), (2, 3), (2, 4), (3, 4), (1, 2), . . . .

Das zyklische Jacobi-Verfahren konvergiert ebenfalls quadratisch ([Wilkinson, 1962],[van Kempen, 1966]). Da es ohne quadratischen Suchaufwand auskommt konvergiertes jedoch deutlich schneller.


Jacobi-VerfahrenParallele Ausführung

Beide TeilschritteA 7→ A′ = JTA, A′ 7→ A′J

eines Jacobi-Schrittes bestehen aus der Aufdatierung von jeweils zwei Zeilen bzw.Spalten. Somit können zwei oder mehr Jacobi-Teilschritte mit disjunkten Index-paaren gleichzeitig vorgenommen werden, da diese jeweils unabhängig voneinandersind.Im o.g. Fall n = 4 könnten die 6 verschiedenen Schritte im zyklischen Jacobi-Verfahren in drei Gruppen aufgeteilt werden:

{(1, 2), (3, 4)}, {(1, 3), (2, 4)}, {(1, 4), (2, 3)}

in denen zwei Prozessoren jeweils zwei Schritte parallel ausführen können. Ist ndeutlich größer als die Anzahl Prozessoren, können auch Block-Verfahren eingesetztwerden.


Jacobi-VerfahrenWeitere Bemerkungen

Trotz der quadratischen Konvergenz ist das Jacobi-Verfahren im Allgemeinenlangsamer als QR.Ist eine approximative Eigenvektor-Matrix V bereits gegeben, so istVTAVbereits fast diagonal, was vom Jacobi-Verfahren (nicht aber vomQR-Verfahren) ausgenutzt werden kann.Für A positiv definit haben [Demmel & Veselić, 1992] gezeigt, dass dasJacobi-Verfahren alle Eigenwerte mit großer relativer Genauigkeitapproximieren kann.


Inhalt

1 Organisatorisches

2 Einleitung

3 Das QR-Verfahren




Bisektionsverfahren

Idee: Eigenwertbestimmung für symmetrische Tridiagonalmatrix durch Nullstellen-suche für charakteristisches Polynom p(λ) = det(A− λI ) auf der reellen Achse.

A =

a1 b1

b1 a2. . .

. . . . . . bn−1bn−1 an

∈ Rn×n, bj 6= 0 ∀j = 1, . . . , n− 1.

Nullstellen eines Polynoms, welches durch deine Koeffizienzen gegeben ist, imAllgemeinen schlecht konditioniert. Koeffizienten der Tridiagonalmatrixbestimmen die Nullstellen (Eigenwerte) jedoch mitweniger Sensitivität.Besonders geeignet, wenn nur Eigenwerte in einem Teil des Spektrums, etwanur größte Eigenwerte oder solche in einem gegebenen reellen Intervall, vonInteresse.Tridiagonalstruktur kann entscheidend ausgenutzt werden.Führt auf Verfahren mit Aufwand O(n) pro Eigenwert (vgl. O(n2) beimQR-Verfahren).


BisektionsverfahrenTrennungseigenschaft

Wesentlich geht bei Bisektionsverfahren die Trennungseigenschaft (engl. interlacingproperty) der Eigenwerte Hermitescher Matrizen ein (vgl. Folie 146).

Bezeichnung: Sei A(k) ∈ Rk×k die k-te führende Hauptuntermatrix von A mit(den notwendig paarweise verschiedenen) Eigenwerten

Λ(A(k)) ={λ(k)1 < λ

(k)2 < · · · < λ

(k)k

}.

Dann gilt

λ(k+1)j < λ

(k)j < λ

(k+1)j+1 , k = 1, . . . , n− 1; j = 1, . . . , k − 1.

Die Trennungseigenschaft macht es möglich, die Anzahl der in einem gegebenenreellen Intervall enthaltenen Eigenwerte von A zu bestimmen.


BisektionsverfahrenTrennungseigenschaft: Beispiel

A =

1 11 0 1

1 2 11 −1

−1.5 −1 −0.5 0 0.5 1 1.5 2 2.5 3

Λ(A( 4))

Λ(A( 3))

Λ(A( 2))

Λ(A( 1))


BisektionsverfahrenTrennungseigenschaft: Beispiel

Aus der Information

detA(1) = 1, detA(2) = −1, detA(3) = −3, detA(4) = 4,

folgt mit derTrennungseigenschaft sukzessiveA(1) besitzt keine negative Eigenwerte.A(2) besitzt genau einen negativen Eigenwert.A(3) besitzt genau einen negativen Eigenwert.A(4) besitzt genau zwei negative Eigenwerte.


BisektionsverfahrenVorzeichenwechsel

Allgemein: bei einer Hermiteschen Tridiagonalmatrix ist die Anzahl negativer Ei-genwerte gegeben durch die Anzahl Vorzeichenwechsel in der Folge

1,detA(1),detA(2), . . . ,detA(n).

Treten in dieser Folge Nullen auf, bleibt diese Aussage richtig, wenn man als Vor-zeichenwechsel die Übergänge

+/0→ −, 0/− → +, aber nicht + /− → 0.

zählt.

Durch einen ShiftA−aI lässt sich so auch die Anzahl Eigenwerte vonA bestimmen,die kleiner als a ∈ R sind.

Die Anzahl Eigenwerte in [a, b) ergibt sich schliesslich aus der Differenz der Anzahlin (−∞, b) und in (−∞, a).


BisektionsverfahrenRekursive Determinantenberechnung

Entwicklung von detA(k) nach der k-ten Zeilen ergibt die rekursive Beziehung

detA(k) = ak detA(k−1) − b2k−1 detA(k−2).

Mit Shift λ ∈ R gilt somit für das charakteristische Polynom

pk(λ) = det(A(k) − λI )

der k-ten Hauptuntermatrix von A− λI die Rekursionsbeziehung

pk(λ) = (ak − λ) pk−1(λ)− b2k−1 pk−2(λ), k = 1, . . . , n,

mit der Initialisierungp−1(λ) ≡ 0, p0(λ) ≡ 1.

Anwendung dieser Iteration für eine Folge von λ-Werten lassen sich durch BisektionEigenwerte in beliebig vorgegebenen Intervallen beliebig genau einschließen.


Inhalt

1 Organisatorisches

2 Einleitung

3 Das QR-Verfahren




Divide & Conquer - VerfahrenGrundidee

Idee: Zerlege symmetrische Tridiagonalmatrix

T =

a1 b1

b1 a2. . .

. . . . . . bn−1bn−1 an

∈ Rn×n

in

T =

[T1 OO T2

]+ bmvvT, m ≈ n

2, v = em + em+1.

Divide-and-Conquer (DQ) Strategie:Löse Eigenproblem für T1, T2 (zwei Probleme der halben Größe).Kombiniere diese Lösungen zur Lösung des Eigenproblems für T .Setze rekursiv fort.


Divide & Conquer - VerfahrenFormulierung mittels Spektralbasis

Spektralkoordinaten: seien

Ti = QiΛiQTi , QT

i Qi = I , i = 1, 2.

Dann ist

T =

[Q1Λ1Q

T1 O

O Q2Λ2QT2

]+ bmvvT

=

[Q1 OO Q2

]([Λ1 OO Λ2

]+ bmuuT

)[QT

1 OO QT

2

], u =

[QT

1 OO QT

2

]v .

Setze

D :=

[Λ1 OO Λ2

]=: diag(d1, d2, . . . , dn), d1 ≥ d2 ≥ · · · ≥ dn

sowie ρ := bm. Dann istΛ(T ) = Λ(D + ρuuT).


Divide & Conquer - VerfahrenSäkulargleichung

Somit gilt λ ∈ Λ(T ) = Λ(D + ρuuT) genau dann, wenn

0 = det(D + ρuuT − λI ) = det((D − λI )(I + ρ(D − λI )−1uuT)

).

Lemma 4.2Für zwei Vektoren x ,y ∈ Rn gilt

det(I + xyT) = 1 + yTx .

Somit istdet(I + ρ(D − λI )−1uuT

)= 1 + ρuT(D − λI )−1u

d.h. die Eigenwerte von T sind die Nullstellen von

f(λ) = 1 + ρ

n∑i=1

u2idi − λ

.

Die Gleichung f(λ) = 0 heißt auch Säkulargleichung (engl. secular equation).


Divide & Conquer - VerfahrenGestalt von f(λ) (Beispiel)

−1 0 1 2 3 4 5 6−6

−4

−2

0

2

4

6

f(λ) = 1 +1

2

(1

1− λ+

1

2− λ+

1

3− λ+

1

4− λ

)Oliver Ernst (Numerische Mathematik) Numerische Lineare Algebra Wintersemester 2013/14 96 / 177

Divide & Conquer - VerfahrenBestimmung der Nullstellen (Grundidee)

Wegen

f ′(λ) = ρ

n∑i=1

u2idi − λ

> 0, λ 6∈ Λ(D),

ist f (für ρ > 0) streng monoton wachsend auf den Intervallen

(−∞, dn), (dn, dn−1), . . . , (d2, d1), (d1,∞).

Die Nullstellen von f werden getrennt durch {di}ni=1. Eine weitere Nullstelleliegt rechts von d1 (ρ > 0) bzw. links von dn (ρ < 0).f ist außerhalb von Λ(D) beliebig oft differenzierbar.Das Newton-Verfahren mit Startwert in jeweils einem Intervall konvergiert ineiner vom Eigenwert unabhängigen Anzahl von Schritten.Da die Auswertung f und f ′ jeweils O(n) Flops erfordert ergibt dies O(n)Flops pro Eigenwert, d.h. O(n2) Flops zur Berechnung aller Eigenwerte vonD + ρuuT.


Divide & Conquer - VerfahrenEigenvektoren

Die Eigenvektoren lassen sich einfach aus den Eigenwerten bestimmen:

Lemma 4.3Ist α Eigenwert von D + ρuuT, so ist (D − αI )−1u ein zugehöriger Eigenvektor.

Da D−αI diagonal erfordert die Berechnung eines solchen Eigenvektors nur O(n)Flops.Diese Formel ist allerdings numerisch instabil – eine stabile Implementierung ist imBuch [Demmel, 1997] nachzulesen.


Divide & Conquer - VerfahrenAlgorithmus

Algorithmus 4 : DQ-Verfahren für das symmetrisch-tridiagonale Eigenpro-blem.Gegeben : T ∈ Rn×n,T = TT tridiagonal.

1 if T ∈ R1×1 then2 return Q = 1, Λ = T

3 else

4 Zerlege T =

[T1 OO T2

]+ bmvvT

5 (Λ1,Q1) aus rekursivem Aufruf für (T1)6 (Λ2,Q2) aus rekursivem Aufruf für (T2)

7 Bilde D + ρuu> aus Λ1,Λ2,Q1,Q2

8 Bestimme Eigenwerte Λ und Eigenvektoren Q ′ von D + ρuuT

9 Berechne Eigenvektoren von T als Q =

[Q1 OO Q2

]Q ′

10 return Q , Λ


Divide & Conquer - VerfahrenKomplexität

Bezeichnet t(n) die zur Anwendung von Algorithmus 4 erfordeliche Anzahl Flops,so ergibt sich

t(n) = 2t(n/2) 2 rekursive Aufrufe für kleinere Eigenprobleme

+O(n2) Eigenwerte von D + ρuuT

+O(n2) Eigenvektoren von D + ρuuT

+ cn3 MultiplikationQ =

[Q1 OO Q2

]Q ′

Ignoriert man die O(n2) Terme ergibt dies t(n) = 2t(n/2) + cn3. Auflösungder Rekursion ergibt t(n) ≈ c 43n

3 (ÜA).Für voll besetzte Matrizen Q1,Q2,Q

′ erhält man mit StandardMatrix-Matrix Multiplikation im letzten Schritt als Konstante c = 1.Durch Deflation erhält man für Q ′ eine dünn-besetzte Matrix.Mit Hilfe der schnellen Multipolmethode kann die Komplexität des Verfahrensauf O(n logp n) reduziert werden (siehe [Demmel, 1997], [Gu & Eisenstat,1995]).


Divide & Conquer - VerfahrenDeflation

Liegt nicht der generische Fall (alle di paarweise verschieden und alle ui 6= 0)vor, so besitzt f nur k < n vertikale Asymptoten und somit k < n Nullstellen.Man kann jedoch zeigen (ÜA), dass die fehlenden n− k Eigenwerte vonD + ρuuT gegeben sind durch {di : ui = 0 oder di = di+1}.Diese Phänomen wird Deflation genannt.In der Praxis wird diese vorgenommen, wenn di+1 ≈ di oder wenn |ui|hinreichend klein.Der belschleunigende Effekt von Deflation besteht darin, dass im Fall ui = 0der Vektor ei zugehöriger Eigenvektor ist. (ÜA)Somit ist Ei auch die i-te Spalte von Q ′ und zur Multiplikation mit Q1 oderQ2 mit dieser Spalte fällt kein Aufwand an.Eine ahnliche Vereinfachung ergibt sich für di = di+1.


Divide & Conquer - VerfahrenLösung der Säkulargleichung

Newton-Verfahren zur Lösung der Säkulargleichung f(λ) = 0 mitStartwert λj :1 Approximiere f in Umgebung von λj durch Tangente.2 Nehme Nullstelle λj+1 der Tangente als nächste Approximation.

Die Abbildungen auf der nächsten Folie zeigen den Verlauf der Funktion f wennu2i = 1/2 (früheres Bild) durch i2i = 103 ersetzt wird, was immer noch zu groß füreine Deflation ist.

Da f in diesem Fall fast immer eine nahezu waagrechte Tangente besitzt würdedas klassische Newton-Verfahren für nahezu jeden Startwert λ0 aus dem aktuellenIntervall springen.

Eine erfolgreiche Alternative besteht darin, f in (di+1, , di) durch rationale Aus-drücke der Form

f(λ) ≈ h(λ) :=c1

di − λ+

c2di+1 − λ

+ c3

zu approximieren, welche f an der aktuellen approximativen Nullstelle λj berühren.Details in [Demmel, 1997].Oliver Ernst (Numerische Mathematik) Numerische Lineare Algebra Wintersemester 2013/14 102 / 177

Divide & Conquer - VerfahrenLösung der Säkulargleichung

−1 0 1 2 3 4 5 6−6

−4

−2

0

2

4

6

1.99 1.995 2 2.005 2.01−6

−4

−2

0

2

4

6

f(λ) = 1 + 10−3

(1

1− λ+

1

2− λ+

1

3− λ+

1

4− λ

)


Inhalt

1 Organisatorisches

2 Einleitung

3 Das QR-Verfahren




Berechnung der Singulärwertzerlegung

Bekanntlich existiert für jede Matrix A ∈ Cm×n die Singulärwertzerlegung (engl.singular value decomposition, SVD)

A = UΣVH

mit unitären Matrizen U ∈ Cm×m, V ∈ Cn×n sowie einer „Diagonalmatrix“

Σ =

[Σr OO O

]∈ Rm×n, Σr = diag(σ1, . . . , σr),

σ1 ≥ · · · ≥ σ2 ≥ σr > 0.

Annahmen:Wir beschränken uns hier auf den reellen Fall A ∈ Rm×n. Dann sind U , Vorthogonal.Wir können m ≥ n annehmen, andernfalls betrachten wir AT anstelle von A,aus deren SVD sich sofort die von A ergibt.


Berechnung der SingulärwertzerlegungRückführung auf Eigenwertaufgabe

Aus der SVD von A ergeben sich sofort die Spektralzerlegungen

ATA = VΣTΣVT, AAT = UΣΣTUT

von ATA bzw. AAT, was etwa folgendes Vorgehen nahelegt:

1 Bilde ATA,2 berechne Spektralzerlegung ATA = VΛVT, Λ = diag(λ1, . . . , λn),3 Σr := diag(

√λ1, . . . ,

√λr) aus nichtnegativen Eigenwerten,

4 Berechne U aus UΣ = AV (z.B. QR-Zerlegung).

In endlich genauer Arithmetik geht jedoch im Allgemeinen bei der Bildung von ATAbzw. AAT Genauigkeit unwiederbringlich verloren. Daher sind Verfahren gefragt,welche direkt mit A arbeiten.


Berechnung der SingulärwertzerlegungRückführung auf Eigenwertaufgabe: Matlab-Beispiel

format long e, rng(2)[U,R] = qr(randn(3)); [V,R] = qr(rand(2));S = diag([1, 1.23456789e-10]);A = U(:,1:2) * S * V’;Lambda = eig(A’*A)Sigma = svd(A)fprintf(’Sqrt(|Lambda|) = %d, %d\n’, sqrt(abs(Lambda)))

Lambda = -2.775557561562891e-171.000000000000000e+00

Sigma = 1.000000000000000e+001.234567736930711e-10

Sqrt(|Lambda|) = 5.268356e-09, 1


Berechnung der SingulärwertzerlegungRückführung auf Eigenwertaufgabe: Sensitivität

Sensitivität von Eigenwerten/Singulärwerten:

|λk(ATA + E)− λk(ATA)| ≤ ‖E‖2|σk(A + E)− σk(A)| ≤ ‖E‖2

Rückwärtsstabiles Verfahren liefert:

σ̃k = σk(A + E), ‖E‖/‖A‖ = O(εMach)

Stabile Berechnung von λk(ATA):

|λ̃k − λk| = O(εMach‖ATA‖) = O(εMach‖A‖2)

Wurzel ziehen

|σ̃k − σk| = O

(|λ̃k − λk|√

λk

)= O

(εMach

‖A‖2

σk

)Fazit: Für kleine Singulärwerte σk � ‖A‖ wird Sensitivität faktisch quadriert.


Berechnung der SingulärwertzerlegungErweiterte Matrix

Für m = n ergibt sich mit A = UΣVT für die erweiterte Matrix

H :=

[O AT

A O

]∈ R2n×2n

die Spektralzerlegung[O AT

A O

] [V VU −U

]=

[V VU −U

] [Σ OO −Σ

]mit der (bis auf den Faktor 1/

√2) orthonormalen Eigenvektormatrix

[V VU −U

].

Beobachtungen:

Die Singulärwerte von A sind die Beträge der Eigenwerte von H .Die Singulärvektoren von A können aus den Eigenvektoren vonH extrahiertwerden.


Berechnung der SingulärwertzerlegungErweiterte Matrix: rechteckiger Fall

Für m > n gilt

VTATAV = ΣTΣ = diag(σ21 , . . . , σ

2n),

UTAATU = ΣΣT = diag(σ21 , . . . , σ

2n, 0, . . . , 0︸︷︷︸

m−n

).

Zerlegt man die Spalten von U in

U =[U1 U2

], U1 ∈ Rm×n,U2 ∈ Rm×m−n,

so ist

Q :=1√2

[V V O

U1 −U1

√2U2

]∈ R(m+n)×(m+n)

orthogonal und es gilt

QT

[O AT

A O

]Q = diag(σ1, . . . , σn,−σ1, . . . ,−σn, 0, . . . , 0︸︷︷︸

m−n

).


Berechnung der SingulärwertzerlegungDas Verfahren von Golub & Kahan

Ähnlich wir beim QR-Verfahren ist es vorteilhaft, bei der Berechnung der SVD dieMatrix zunächst (orthogonal) zu transformieren, in diesem Fall in Bidiagonalform.Die Rechnung gliedert sich dann in zwei Phasen:× × × ×× × × ×× × × ×× × × ×× × × ×× × × ×

Phase 1−−−−−−−−−−→

Bidiagonalisierung

× ×× ×× ××

Phase 2−−−−−→Iteration

××××

Phase 1 wird Golub-Kahan-Bidiagonalisierung (nach [Golub & Kahan, 1965])genannt. Diese erfordert mit Householder-Spiegelungen O(mn2) Flops.In Phase 2 wird eines der Eigenwertverfahren für symmetrische Matrizenangewandt, wobei nur mit der Bidiagonalmatrix gearbeitet wird. DasQR-Verfahren benötigt typischerweise O(n log | log εMach|) Iterationsschritte,um jeden Eigenwert bis auf Machinengenauigkeit zu berechnen, bei O(n)Flops pro Schritt also insgesamt O(n2) Iterationen.


Berechnung der SingulärwertzerlegungBidiagonalisierung

Abwechselnde Householder-Reflektionen von links (U>j ·) bzw. rechts (·Vj) führenschrittweise auf Bidiagonalform:

× × × ×× × × ×× × × ×× × × ×× × × ×× × × ×

U>

1 ·−−−→

× × × ×0 × × ×0 × × ×0 × × ×0 × × ×0 × × ×

·V1−−→

× × 0 0× × ×× × ×× × ×× × ×× × ×

U>2 ·−−−→

× ×× × ×0 × ×0 × ×0 × ×0 × ×

·V2−−→

× ×× × 0× ×× ×× ×× ×

U>

3 ·−−−→

× ×× ×× ×0 ×0 ×0 ×

. . .



U>4 ·−−−→

× ×× ×× ××00

=: B =: Q>L AQR

QL = U1U2U3U4,

QR = V1V2

A = QLBQTR

Bidiagonalform ist spätestens nach m− 2 Housholder-Spiegelungen von links sowien− 2 von rechts erreicht.

Nach erfolgter Berechnung der SVD B = UBΣVTB der Bidiagonalmatrix B erhal-

ten wir mitA = QLUB︸︷︷︸

=:U

ΣVTBQ

TR︸︷︷︸

=:VT

die SVD von A.



Mit B = QTLAQR gilt ebenso[

QTR O

O QTL

] [O AT

A O

] [QR OO QL

]=

[O QT

RATQL

QTLAQR O

]=

[O BT

B O

]=: HB ,

d.h. der Golub-Kahan-Bidiagonalisierung entspricht einer orthogonalen Ähnlichkeit-stransformation der Matrix H auf Bandstruktur HB .

Führt man noch eine umordnung der Zeilen und Spalten gemäß der „perfect shuffle“Permutation

1, n+ 1, 2, n+ 2, 3, n+ 3, . . .

durch, so wird HB zu einer symmetrischen Tridiagonalmatrix mit HauptdiagonaleNull und unterer/oberer Nebendiagonale α1, β1, α2, β2, α3 . . . .


Berechnung der SingulärwertzerlegungQR für SVD

Bezeichnung:

B =

α1 β1

α2. . .. . . βn−1

αn

Vereinfachende Annahmen:

B ist quadratisch, andernfalls müssen nur Nullzeilen und -spalten in HB

gestrichen werden.Alle αj sind ungleich Null, andernfalls ist Null ein Singulärwert und kannsofort abgespaltet werden.Alle βj sind ungleich Null, andernfalls kann das Problem in zwei kleinereProbleme zerlegt werden.


Berechnung der SingulärwertzerlegungQR für SVD

Die Matrix HB entspricht der zyklischen Matrix C imProdukt-QR-Eigenwertproblem im Spezialfall k = 2 sowie A1 = AT

2.Daher: QR-Schritte vom Grad 2m mit Shift-Paaren ±µ (ω = −1).Da die beiden Blöcke transponiert zueinander sind, entsprechen sich diebulge-chasing Schritte und können somit an nur einer der beidendurchgeführt werden.m = 1: Golub-Kahan QR-Verfahren für SVD-Berechnung. EntsprichtDoppelshift QR-Verfahren angewandt auf PTHBP (perfect shufflePermutation) mit Shift-Paaren ±µ ohne redundante Operationen.Es gilt: βn−1 → 0 kubisch, αn konvergiert gegen einen Singulärwert von A,der daraufhin deflationiert werden kann.


Berechnung der SingulärwertzerlegungQR für SVD, hohe relative Genauigkeit

[Demmel & Kahan, 1990]: wichtige Verbesserung des Golub-Kahan Verfahrens.Einträge der Bidiagonalmatrix B bestimmen deren Singulärwerte zu hoherGenauigkeit, d.h. winzige relative Störungen der Einträge führen nur zuwinzigen relativen Störungen der Singulärwerte.Dies ist nicht der Fall für die Einträge der Tridiagonalmatrix T = BTB .Idee: Subtraktionsfreie Formulierung des Golub-Kahan QR Schrittes ohneShift (µ = 0). QR-Iteration ohne Shifts bis kleine Singulärwerte gefunden,danach übrige mit geschifteten QR-Schritten.Weitere Verbesserungen mit hoher relativer Genauigkeit für alle Singulärwertesowie Beschleunigung: differential quotient-difference method with shifts(dqds), [Fernando & Parlett, 1994], [Demmel, 1997].


Numerische Lineare Algebra - tu-chemnitz.de · Jacobi-Verfahren Konvergenz,zyklischeIndexwahl...

Documents

Transcript of Numerische Lineare Algebra - tu-chemnitz.de · Jacobi-Verfahren Konvergenz,zyklischeIndexwahl...