7. iterative l¨osung linearer...
Post on 12-Aug-2019
214 Views
Preview:
TRANSCRIPT
7. Iterative Losung
linearer Gleichungssysteme
1
Grundlagen und Wiederholung (1)
Die Grundlagen decken sich mit dem Stoff, der einen Teil des Kapitels
2 - Numerik ausmacht und bereits in Mathematik behandelt wurde.
Eine zentrale Rolle bei numerischen Berechnungen spielen lineare Glei-
chungssysteme
• Es sind die am haufigsten auftretenden numerischen Probleme
• Anwendungsgebiete sind z.B.
* fast alle naturwissenschaftlich-technischen Problemstellungen
vom Wetterbericht bis zur Warmeentwicklung auf einer Koch-
platte oder der Planung der Leiterbahnen auf Mikrochips,
* Bildverarbeitung oder z.B Beleuchtungsprobleme in der Com-
putergrafik,
* wirtschaftlichen Fragestellungen wie Versicherungskosten oder
Borsenkursvorhersage.
2
Grundlagen und Wiederholung (2)
Losungsverfahren
• Die direkten Verfahren liefern eine mit Rundungsfehlern behaftete
Losung nach endlich vielen Schritten.
• Die iterativen Verfahren beginnen mit einer Anfangsnaherung und
produzieren eine verbesserte Naherungslosung nach endlich vielen
Schritten.
• Falls moglich wird das Problem mit einem direkten Verfahren be-
rechnet und anschließend werden die Rundungsfehler mit einem
iterativen Verfahren verringert.
3
Grundlagen und Wiederholung (3)
Problemstellung: Berechne den Vektor x = (x1, x2, . . . xn) aus
a1,1x1 + a1,2x2 + · · · a1,nxn = b1
a2,1x1 + a2,2x2 + · · · a2,nxn = b2
·
·
an,1x1 + an,2x2 + · · · an,nxn = bn
oder in Matrix-Schreibweise
Ax = b
Bemerkung: Vektoren werden hier ohne Vektorpfeil geschrieben
4
Wiederholung Gauß-Verfahren (1)
Schulbeispiel:
E1 : x1 + x2 + 3x4 = 4E2 : − x2 − x3 − 5x4 = −7E3 : − 4x2 − x3 − 7x4 = −15E4 : 3x2 + 3x3 + 2x4 = 8
Mit Matrix:
A =
1 1 0 3
0 −1 −1 −5
0 −4 −1 −7
0 3 3 2
, x =
x1
x2
x3
x4
, b =
4
−7
−15
8
5
Wiederholung Gauß-Verfahren (2)
Erlaubte Transformationen zur Losung des Gleichungssystems
• Multiplizieren einer Zeile (Gleichung) mit einer Zahl verschieden
von Null
• Addieren eines Vielfachen einer Zeile zu einer anderen Zeile
• Vertauschen von Zeilen (Gleichungen) bzw. Spalten (Unbekann-
ten, entspricht Umnummerierung)
Mit Hilfe dieser Transformationen reduziere Gleichungssystem auf ein
Dreieckssystem.
6
Wiederholung Gauß-Verfahren (3)
Zweite Spalte des Schulbeispiels:
1 1 0 3
0 −1 −1 −5
0 −4 −1 −7
0 3 3 2
·
x1
x2
x3
x4
=
4
−7
−15
8
·
a3,i = a3,i − a3,2/a2,2 · a2,i b3 = b3 − a3,2/a2,2 · b2
a4,i = a4,i − a4,2/a2,2 · a2,i b4 = b4 − a4,2/a2,2 · b2
7
Wiederholung Gauß-Verfahren (4)
Schulbeispiel wird zu
E1 : x1 + x2 + 3x4 = 4E2 : − x2 − x3 − 5x4 = −7E3 : 3x3 + 13x4 = 13E4 : − 13x4 = −13
oder allgemein
a1,1 a1,2 · · · a1,n
a2,2 · · · a2,n
. . . ...
an,n
·
x1
x2
...
xn
=
b1
b2
...
bn
·
8
Wiederholung Gauß-Verfahren (5)
Ein Dreieckssystem ist leicht zu losen. Aus E4:
x4 = 1
x4 einsetzten in E3:
3x3 +13 = 13 → x3 = 0
x3, x4 einsetzten in E2:
−x2 − 5 = −7 → x2 = 2
x2, x3, x4 einsetzten in E1:
x1 +2+ 3 = 4 → x1 = −1
Der Gauß-Algorithmus ist von der Ordnung O(n3).
9
Norm von Vektoren und Matrizen
Unterschiedliche Definitionen von Normen sind moglich, hier nur die
sogenannte 2-Norm:
• 2-Norm oder euklidische Norm von Vektoren (Lange eines Vek-
tors):
||x||2 =
√
√
√
√
n∑
i=1
x2i
• 2-Norm oder Spektralnorm von Matrizen (Wurzel aus der Summe
der Quadrate der Diagonalelemente bei einer Diagonalmatrix):
||A||2 = maxx6=0||Ax||
||x||=
√
ρ(ATA)
10
Iterative Verfahren (1)
Da das Gauß-Verfahren
• O(n3) ist und damit oft zu lange dauert,
• viele Matrizen nur dunn besetzt sind, das Gauß-Verfahren aber die
Matrix fullt
• und das Gauß-Verfahren zu viele Rundungsfehler hat,
werden meist iterative Gleichungsloser verwendet. Beispiele:
• alle Programme fur Computergraphiken, wie z.B. PovRay,
• alle Programme fur Computersimulationen wie im IMH.
11
Iterative Verfahren (2)
Idee: Fixpunktgleichung
Eine Gleichung der Form F(x) = x heißt Fixpunktgleichung. Ihre
Losungen, also der Werte, der die Gleichung F(x) = x erfullen, heißen
Fixpunkte
Definition:
Gegeben sei F : [a, b] → R,
x0 ∈ [a, b]. Die rekursive Folge
xn+1 := F(xn), n = 0,1, . . .
heißt Fixpunktiteration von F
zum Startwert x0.
X0
X1
X2
X3
X1
X2
X3
12
Iterative Verfahren (3)
Anwendungsbeispiel
Bestimme die Nullstelle von p(x) = x3−x+0.3 oder lose die Gleichung
x = x3 +0.3.
Methode: Fuhre die Fixpunktiteration xn+1 = F(xn) = x3n+0.3, durch
Ergebnis:
• Startwert x0 = 0 : konvergiert gegen x = 0.3389 . . .
• Startwert x0 = 1: divergiert
Z.B. kann die Gleichung x = 2sin(x) + 1 nur iterativ gelost werden!
13
Iterative Verfahren (4)
Iterative Losungsverfahren konstruieren ganz allgemein eine Losung
aus einer Startlosung:
x(0) → x(1) → x(2) → . . .
Angewendet auf die Berechnung der Losung eines linearen Gleichungs-
systems ist die Konvergenz abhangig von den Eigenschaften der Matrix
und den Details des iterativen Losungsverfahrens.
Es gibt zwei Gruppen von Verfahren
• die “klassischen” Verfahren, das Jacobi- und das Gauß-Seidel-
Iterationsverfahren und darauf aufbauend Relaxationsverfah-
ren, entwickelt im spaten 18. Jahrhundert, werden heute noch
angewandt.
• Krylov Unterraum-Methoden, die allgemeiner und oft schneller
sind, jedoch ist in vielen Fallen die Konvergenz unklar.
14
Richardson Iterationsverfahren (1)
Richardson Iterationsverfahren:
Idee: Formuliere ein passendes Fixpunktproblem
b = Ax = (A− I)x+ x ⇔ x = b+ (I− A)x := F(x)
Iterationsvorschrift:
x(t+1) = b+ (I− A)x(t) = x(t) + (b−Ax(t)) = x(t) + r(t)
Fur x(t+1) = x(t) gilt r(t) = 0 oder b− Ax(t) = 0
Der gesuchte Vektor x ist der Fixpunkt der Gleichung und wir oft auch
als x bezeichnet.
Der Vektor r(t) heißt Rest- oder Residuenvektor (manchmal Residu-
umsvektor oder kurz Residuum) und ist ein Maß fur den Fehler.
15
Richardson Iterationsverfahren (2)
Abbruchbedingung der Iterationsverfahren: Der Betrag des Re-
siduenvektors ||r(t)|| = ||b − Ax(t)|| = ||A(x − x(t))|| oder besser der
relative Betrag ||r(t)||/||x(t)|| muss klein sein.
||x|| steht fur die 2-Norm oder der Lange des Vektors x
||x|| =
√
√
√
√
n∑
i=1
x2i
Fehlerbetrachtung im t-ten Schritt
x− x(t) = x− x(t−1) − (b− Ax(t−1))
= x− x(t−1) − (Ax− Ax(t−1))
= (I− A)(x− x(t−1))
16
Richardson Iterationsverfahren (3)
Mit der 2-Norm
||x− x(t)||2 = ||(I− A)(x− x(t−1))||2
≤ ||(I− A)||2 · ||(x− x(t−1))||2
≤ ||(I− A)||22 · ||(x− x(t−2))||2
· · ·
≤ ||(I− A)||t2 · ||(x− x(0))||2
Konvergenz liegt sicher vor fur
||(I− A)|| < 1
oder A ist fast eine Einheitsmatrix.
(Die 2-Norm der Matrix ist die “Spektralnorm” und wird hier nicht
naher betrachtet, siehe Kapitel 2 zu Matrixnorm)
17
Jacobi Iterationsverfahren (1)
Verbesserung:
In vielen praktischen Problemen sind die Matrixeintrage auf der Dia-
gonalen groß. Betrachte dann das modifizierte Problem
D−1Ax = D−1b mit D = diag(A)
Diese Anderung andert die Losung nicht!
D−1A ≈ I bzw. ||(I−D−1A)|| < 1 ist aber eher erfullt als ||(I−A)|| < 1
Beispiel:
A =
5 2 12 7 4−1 2 8
, D−1 =
1/5 0 00 1/7 00 0 1/8
, I−D−1A =
0 2/7 1/82/5 0 4/8−1/5 2/7 0
,
Matrizen, deren großte Eintrage auf der Diagonalen sind, kommen z.B.
bei Beleuchtungsproblemen (Kapitel 2) oder bei sogenannten partiel-
len Differentialgleichungen vor (Kapitel 10).
18
Jacobi Iterationsverfahren (2)
Aus Ax = b wurde D−1Ax = D−1b, also, ersetze im Richardson-
Verfahren A → D−1A und b → D−1b:
x(t+1) = x(t) + (D−1b−D−1Ax(t)) oder
Dx(t+1) = Dx(t) + (b−Ax(t)) = b+ (D − A)x(t)
Das ist das Jacobiverfahren
Historisch wurde das Verfahren anders hergeleitet, so wie es auch meist
in der Literatur zu finden ist.
Ersetze Ax durch (D + A−D)x und forme den Ausdruck in eine Fix-
punktgleichung um
Ax = Dx+ (A−D)x = b
Das fuhrt zur selben Iterationsvorschrift:
Dx(t+1) = b+ (D −A)x(t)
19
Jacobi Iterationsverfahren (3)
Elementweise geschrieben:
ajjx(t+1)j = bj + ajjx
(t)j −
∑
k
ajkx(t)k
x(t+1)j =
1
ajj
bj −∑
k 6=j
ajkx(t)k .
• Multipliziere die Matrix ohne Diagonalelemente mit x(t)
• Ziehe das Ergebnis vom Vektor b ab
• Teile jedes Element durch das entsprechende Diagonalelement
• Fuhre das solange durch, bis der Residuenvektor klein ist
20
Gauß-Seidel Verfahren (1)
Weitere Verbesserung: Zerlege die Matrix in 3 Teile:
• D: Diagonalteil von A
• −L: Unterdiagonalteil von A
• −U : Oberdiagonalteil von A
Lose
b = Ax = (D − L− U)x = (D − L)x− Ux
uber die Iterationsvorschrift
(D − L)x(t+1) = b+ Ux(t) = b− (A− (D − L))x(t)
= (D − L)x(t) + (b−Ax(t)) = (D − L)x(t) + r(t)
21
Gauß-Seidel Verfahren (2)
Das Gauß-Seidel Verfahren lautet
x(t+1) = x(t) + (D − L)−1(b− Ax(t))
und ist gleich einer Richardson Iteration angewandt auf
(D − L)−1Ax = (D − L)−1b
Das Verfahren ist konvergent, falls ||I− (D − L)−1A||2 < 1
In vielen Anwendungen, basierend auf der diskretisierten Form soge-
nannter partieller Differentialgleichungen, kommt eine Variante dieser
Methode zum Einsatz (Kapitel 10).
22
Gauß-Seidel Verfahren (3)
Elementweise geschrieben:
(D − L)x(t+1) = b+ Ux(t)
ajjx(t+1)j +
∑
k<j
ajkx(t+1)k = bj −
∑
k>j
ajkx(t)k .
Berechnet x(t+1)1 , und damit x
(t+1)2 usw.
a11x(t+1)1 = b1 −
∑
k>1
a1kx(t)k
a22x(t+1)2 + a21x
(t+1)1 = b2 −
∑
k>2
a2kx(t)k
a33x(t+1)3 + a31x
(t+1)1 + a32x
(t+1)2 = b3 −
∑
k>3
a3kx(t)k
... = ...
23
Relaxationsverfahren
Die Wahl der Matrix, mit der die ursprunglich Fixpunktgleichung mo-
difiziert wird, ist im Prinzip beliebig (siehe Prakonditionierung)
Das SOR (successive over-relaxation) Verfahren modifiziert das Gauß-
Seidel Verfahren durch
(D − L) → (D
ω− L)
Eingesetzt und multipliziert mit ω folgt
(D − ωL) x(t+1) = [(1− ω)D + ωU ]x(t) + ωb
Die richtige Wahl vom ω kann den Algorithmus stark beschleunigen.
24
Prakonditionierung
Verallgemeinerung der obigen Methoden: Anstatt das Originalproblem
zu losen, betrachte
M−1Ax = M−1b
M sollte leicht zu invertieren sein und in der Nahe von A liegen, damit
M−1A leicht zu berechnen ist und in der Nahe der Einheitsmatrix liegt.
Die Richardson Iterationsvorschrift fur das prakonditionierte System
lautet
x(t+1) = (1−M−1A)x(t)+M−1b = x(t)+M−1(b−Ax(t)) = x(t)+M−1r(t)
mit hoffentlich kleiner Iterationsmatrix C = 1 − M−1A. Jacobi bzw.
Gauß-Seidel Verfahren entsprechen den Prakonditionierungsmatrizen
M = D bzw. M = D − L.
25
Krylov Unterraum-Methoden (1)
Problem:
Die bis jetzt vorgestellten Methoden konvergieren nur, wenn die Itera-
tionsmatrix ||C|| = ||1−M−1A|| klein ist.
Fur das Jacobi-, Gauß-Seidel- und SOR-Verfahren bedeutet das, dass
die Nebendiagonalelemente von A im Vergleich zu den Diagonalele-
menten klein sein mussen.
Krylov Unterraum-Methoden beruhen auf einer Verallgemeinerung
der bisher vorgestellten Fixpunktgleichungen. Aus dem prakonditio-
niertem Richardson Iterationsverfahren wird
x(t+1) = x(t) +M−1r(t) ⇒ x(t+1) = x(t) + α(t)p(t).
Je nach Wahl von α(t) und p(t) ergeben sich unterschiedliche Metho-
den, wobei nur sicher gestellt werden muss, dass x(t) jede mogliche
Richtung annehmen kann, so dass die Losung gefunden werden kann,
und dass der Betrag des Residuums ||b− Axk|| immer kleiner wird.
26
Krylov Unterraum-Methoden (2)
Was ist ein Krylov Unterraum?
Beispiel: Fur die einfachste Iterationsvorschrift: Die Richardson Itera-
tion
x(t+1) = x(t) + r(t).
gilt
r(t) = b−Ax(t)
= b−A(x(t−1) + r(t−1))
= b−A(x(t−1) + b−Ax(t−1))
= (I−A)(b− Ax(t−1))
= (I−A)r(t−1)
= (I−A)2r(t−2)
= . . . = (I− A)tr(0)
27
Krylov Unterraum-Methoden (3)
Mit dem Startvektor x(0) = 0 und damit r(0) = b−Ax0 = b folgt
x(t+1) = x(t) + r(t)
= x(t−1) +t∑
j=t−1
r(j)
= x(0) +t∑
j=0
r(j)
=t∑
j=0
(I−A)j b
28
Krylov Unterraum-Methoden (4)
Die Iterationslosung x(t) =∑t−1
j=1(I− A)j b ist somit Element des Un-
terraums
x(t) ∈ {b, Ab, . . . , At−1b} = Kt(A,b)
Dieser Raum wird als Krylov Unterraum der Dimension t, Kt(A,b)
bezeichnet.
Die Iterationslosung xt liegt bei einem Krylov-Unterraumverfahren in
Kt(A,b).
Je nachdem, wie die Iterationslosung innerhalb des Unterraums
bestimmt wird, gibt es verschiedene Iterations-Strategien.
29
Krylov Unterraum-Methoden (5)
Die bekanntesten sind
• Ritz-Galerkin Ansatz: r(t) orthogonal zu Kt(A,b) (z.B. CG).
Das bedeutet genauer:
xt ∈ Kt und rt = (b−Axt) ∈ K⊥t
• Petrov-Galerkin Ansatz: r(t) orthogonal zu einem allgemeinen Un-
terraum Lt (z.B. BI-CG oder QMR))
• Minimierung von ||r(t)||2 (z.B. MINRES oder GMRES)
• Hybride Verfahren (z.B. Bi-CGSTAB)
• Minimierung des Fehlers ||x− x(t)||2 (z.B. GMERR)
Es werden laufend neue Verfahren, angepasst an spezielle Probleme
entwickelt.
30
Einschub: Schreibweisen fur Vektoren und Matrizen
• Vektoren werden mit einem Vektorpfeil nur geschrieben, wenn es
sich um Vektoren im Raum handelt, ansonsten meist durch “fetten
Schriftsatz” gekennzeichnet.
• Ein Skalarprodukt zweier Vektoren wird geschrieben als
* a◦b oder a·b, haufig bei Ingenieuren und in der Schule verwendet
* 〈a,b〉 oder (a,b), Bra-Ket Schreibweise, beliebt bei z.B. Physi-
kern und angewandten Mathematikern
* aTb, Matrix-Schreibweise, beliebt bei Mathematikern
• Dementsprechend gilt z.B.
a ◦ (Mb) ≡ aTMb ≡ (a,Mb)
31
CG (1)
Das konjugierte Gradientenverfahren (CG) ist Grundlage vieler moder-
ner Iterationsverfahren (entwickelt 1952 von Stiefel und Hestens).
• Es konvergiert (falls keine Rundungsfehler vorliegen) fur positiv
definite und symmetrische Matrizen
a(x) = xTAx = (x, Ax) > 0; Ai,j = Aj,i
der Große n × n in n Schritten und ist somit eine schnelle und
sichere Methode zur Losung des Gleichungssystems.
• Die Rundungsfehler fuhren bei großen Systemen jedoch dazu, dass
man in vielen Fallen nicht n sondern ca. 3n Schritte anwendet.
• Es ist das einzige iterative Verfahren, fur das die Konvergenz im
Allgemeinen bewiesen wurde.
• Da jeder Schritt einer Matrixmultiplikation entspricht, lohnt sich
das Verfahren besonders fur dunn besetzte Matrizen.
32
CG (2)
Grundlage ist ein Optimierungsproblem. Die Losung von
Ax = b
ist ein Minimum der Funktion
f(x) =1
2(x, Ax)− (b,x)
und umgekehrt. Begrundung: Sie x der Losungsvektor von Ax = b, so
gilt fur einen beliebigen Vektor x+ p
f(x+ p) =1
2(x+ p, A(x+ p))− (b,x+ p)
= f(x) + (p, (Ax− b)) +1
2(p, Ap) = f(x) +
1
2(p, Ap) > f(x)
Das Verfahren setzt sich aus 2 Teilen zusammen: die Methode des
“steilsten Abstiegs” und die Methode der “konjugierten Richtung”.
33
CG (3)
Die Methode des “steepest descent” der Gradientenverfahren:
• Richtung: Beginnt man mit einem Vektor x+ p und mochte ent-lang des steilsten Abstiegs das Minimum erreichen, so muss manentlang des Negativen der Ableitung der Funktion f(x) nach x
gehen, also entlang
−f′(x) = −(Ax− b) = r.
Das Residuum gibt die Richtung des steilsten Abstiegs an.
• Relaxationsfaktor: Betrachte die Funktion f(x+αp) und bestimmtfur einen Vektor p das Minimum im Parameter α:
df(x+ αp)
dα=
d
dα(f(x) + (αp, (Ax− b)) +
1
2α2(p, Ap))
= (p, (Ax− b)) + α(p, Ap) = 0
oder
αopt =(p, r)
(p, Ap)
34
CG (4)
Algorithmus: Ausgangspunkt ist eine verallgemeinerte Richardson Ite-
rationsvorschrift x(i+1) = x(i) + αip(i).
• Verwende die Richtung des steilsten Abstiegs p = r und αopt.
• Setze x(0) und berechne r(0) = b−Ax(0)
• Fuhre folgende Schritte durch:
αi =(r(i), r(i))
(r(i), Ar(i))
x(i+1) = x(i) + αir(i)
r(i+1) = b−Ax(i+1) = r(i) − αiAr(i)
Das Verfahren konvergiert immer fur positiv definite und symmetri-
sche Matrizen, aber meist sehr langsam, da die Richtung des steilsten
Abstiegs zu einer oszillierenden Bewegung der Losung fuhren kann.
35
CG (5)
Beispiel aus G.Opfer:
Bestimme die Losung des Gleichungssystem(
1 00 10
)
·
(
x1x2
)
=
(
00
)
mit der Methode des steilsten Abstiegs. Das aquivalente System ist:
Bestimme das Minimum der Funktion
f(x) =1
2(x, Ax)− (b,x)
=1
2
(
x1 x2)
·
(
1 00 10
)
·
(
x1x2
)
−(
0 0)
·
(
x1x2
)
=1
2(x21 +10x22)
36
CG (6)
Wahle als Startvektor
x(0)T= (1.0,1.0)
Der Losungsvektor
oszilliert zur Losung
-0.1
-0.08
-0.06
-0.04
-0.02
0
0.02
0.04
0.06
0.08
0.1
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
min(x_1^2 + 10 x_2^2)
37
CG (7)
Die Methode der “konjugierten Richtung”:
Ausgangspunkt fur Methode des “steepest descent” war die Funktion
f . Als Richtung wurde der Gradient festgelegt und die Lange des
Gradientenvektors durch das Minimum der Funktion f
f(x+ αp) = f(x) + α(p, (Ax− b)) +1
2α2 (p, Ap)
= f(x) + α(p, Ax) +1
2α2 (p, Ap)− α(p,b)
als Funktion von α bestimmt. Verfahren zur Vermeidung der Oszilla-
tionen (nur zum Verstandnis, ohne Beweise):
• Die Richtungen, in der sich x(i) entwickelt, seien in einem Unter-
raum P (i) =< p(0), p(1), . . . p(i−1) >
38
CG (8)
• Der Vektor x(i) lasst sich dann schreiben als
x(i) = x(i−1) + α(i−1)p(i−1) = x(0) +i−1∑
j=0
α(j)p(j)
• Wahle als nachste Richtung zur Berechnung von x(i+1) nicht den
Gradienten, sondern eine Richtung p(i), so dass der 2. Summand
in f(x+ αp) fur alle Betrage von x in Richtung x(i) wegfallt
(p(i), Ax(i)) = 0
• Da x(i) eine Linearkombination der x(j), j = 1, . . . , i−1 ist, bedeutet
das
(p(i), Ap(j)) = 0 fur j ≤ i
39
CG (9)
• Die Vektoren p(0),p(1), . . . ,p(i) heißen paarweise zu A-konjugierte
Richtungen, da gilt
(p(k), Ap(j)) = 0 fur k 6= j
• Bestimme nun, wie gehabt, das optimale α.
αopt =(p(i), r(i))
(p(i), Ap(i))
• und berechne damit
x(i+1) = x(i) + αp(i)
40
CG (10)
Algorithmus: (bis jetzt)
Verwende einen Satz von konjugierten Suchrichtungen pi.
Setze x(0) und r(0) = b− Ax(0)
Dann fuhre folgende Schritte durch:
αi =(r(i),p(i))
(p(i), Ap(i))
x(i+1) = x(i) + αip(i)
r(i+1) = b−Ax(i+1) = r(i) − αiAp(i)
Die verbleibende Aufgabe ist es, die Suchrichtungen p(i) gunstig zu
wahlen, dass sich schnell eine gute Naherung ergibt.
41
CG (11)
Die “konjugierte Gradienten-Methode” von Hestens und Stiefel
Hier werden die konjugierten Richtungen aus den Restvektoren kon-
struiert werden.
p(0) = r(0)
p(i) = r(i) + β(i)p(i−1).
Ist β(i) = 0, ergibt sich die Methode des Gradientenverfahren. Die
Koeffizienten β(i) werden so gewahlt, dass die Richtungen p(i) wie
gefordert zueinander konjugiert sind.
0 = (p(m), Ap(i)) = (r(m), Ap(i))+β(m)(p(m−1), Ap(i)) fur i = 0, . . . ,m−1
42
CG (12)
Da
(p(m−1), Ap(i)) = 0 fur i < m− 1
folgt
β(m−1) = −(r(m), Ap(m−1))
(p(m−1), Ap(m−1))
Ohne Beweis: Jetzt gilt fur den neuen Residuenvektor
(r(i+1),p(j)) = 0 fur j = 0, . . . , i,
d.h. der Vektor steht senkrecht auf dem Raum P (i). Da sich der Vek-
torraum bei jedem Iterationsschritt um eine Dimension vergroßert, ist
sicher gestellt, dass das Verfahren nach n Schritten eine Losung findet.
Nach einigen Umrechnungen ergibt sich CG-Algorithmus.
43
CG (13)
Endgultiger Algorithmus (ohne Beweis!):
Starte mit p(0) = r(0) = b−Ax(0). Dann fuhre folgende Schritte durch:
α(i−1) =|r(i−1)|2
(p(i−1), Ap(i−1))
x(i) = x(i−1) + αi−1p(i−1)
r(i) = r(i−1) − αi−1Ap(i−1)
β(i) =|r(i)|2
|r(i−1)|2
p(i) = r(i) + β(i)pi−1
44
CG (14)
Nochmal das Beispiel
aus G.Opfer(
1 00 10
)
·
(
x1x2
)
=
(
00
)
Wahle als Startvektor
wieder x(0)T= (1.0,1.0)
Der Losungsvektor
oszilliert nicht mehr und
die Losung wird in 2
Schritten erreicht.
-0.1
-0.08
-0.06
-0.04
-0.02
0
0.02
0.04
0.06
0.08
0.1
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
min(x_1^2 + 10 x_2^2)
45
CG (15)
Rechenaufwand pro Iteration
1 Matrixmultiplikation A,
2 Skalarprodukte,
3 AXPYs
Es werden (theoretisch) n Iterationen benotigt.
Einfache Erweiterung fur nicht-symmetrische Matrizen
Ist die Matrix nicht symmetrisch und positiv definite, betrachte
ATAx = AT b
und konstruiere die Losung in Ki(ATA,AT b) unter Verdoppelung der
Anzahl der Matrix-Vektor Multiplikationen (geht besser durch Erwei-
terung des Krylov-Unterraums).
46
top related