7. Iterative L¨ osung
linearer Gleichungssysteme
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 h¨aufigsten auftretenden numerischen Probleme
• Anwendungsgebiete sind z.B.
* fast alle naturwissenschaftlich-technischen Problemstellungen vom Wetterbericht bis zur W¨armeentwicklung 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 B¨orsenkursvorhersage.
Grundlagen und Wiederholung (2)
L¨osungsverfahren
• Die direkten Verfahren liefern eine mit Rundungsfehlern behaftete L¨osung nach endlich vielen Schritten.
• Die iterativen Verfahren beginnen mit einer Anfangsn¨aherung und produzieren eine verbesserte N¨aherungsl¨osung nach endlich vielen Schritten.
• Falls m¨oglich wird das Problem mit einem direkten Verfahren be- rechnet und anschließend werden die Rundungsfehler mit einem iterativen Verfahren verringert.
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
Wiederholung Gauß-Verfahren (1)
Schulbeispiel:
E1 : x1 + x2 + 3x4 = 4 E2 : − x2 − x3 − 5x4 = −7 E3 : − 4x2 − x3 − 7x4 = −15 E4 : 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
Wiederholung Gauß-Verfahren (2)
Erlaubte Transformationen zur L¨osung 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.
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
Wiederholung Gauß-Verfahren (4)
Schulbeispiel wird zu
E1 : x1 + x2 + 3x4 = 4 E2 : − x2 − x3 − 5x4 = −7
E3 : 3x3 + 13x4 = 13
E4 : − 13x4 = −13
oder allgemein
a1,1 a1,2 · · · a1,n a2,2 · · · a2,n
. . . ...
an,n
·
x1 x2 ...
xn
=
b1 b2 ...
bn
·
Wiederholung Gauß-Verfahren (5)
Ein Dreieckssystem ist leicht zu l¨osen. 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).
Norm von Vektoren und Matrizen
Unterschiedliche Definitionen von Normen sind m¨oglich, hier nur die sogenannte 2-Norm:
• 2-Norm oder euklidische Norm von Vektoren (L¨ange eines Vek- tors):
||x||2 =
v u u t
n X 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|| =
q
ρ(ATA)
Iterative Verfahren (1)
Da das Gauß-Verfahren
• O(n3) ist und damit oft zu lange dauert,
• viele Matrizen nur d¨unn besetzt sind, das Gauß-Verfahren aber die Matrix f¨ullt
• und das Gauß-Verfahren zu viele Rundungsfehler hat,
werden meist iterative Gleichungsl¨oser verwendet. Beispiele:
• alle Programme f¨ur Computergraphiken, wie z.B. PovRay,
• alle Programme f¨ur Computersimulationen wie im IMH.
Iterative Verfahren (2)
Idee: Fixpunktgleichung
Eine Gleichung der Form F(x) = x heißt Fixpunktgleichung. Ihre L¨osungen, also der Werte, der die Gleichung F(x) = x erf¨ullen, 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 X
1 X
2 X
3
X1
X2
X3
Iterative Verfahren (3)
Anwendungsbeispiel
Bestimme die Nullstelle von p(x) = x3−x+ 0.3 oder l¨ose die Gleichung x = x3 + 0.3.
Methode: F¨uhre 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 = 2 sin(x) + 1 nur iterativ gel¨ost werden!
Iterative Verfahren (4)
Iterative L¨osungsverfahren konstruieren ganz allgemein eine L¨osung aus einer Startl¨osung:
x(0) → x(1) → x(2) → . . .
Angewendet auf die Berechnung der L¨osung eines linearen Gleichungs- systems ist die Konvergenz abh¨angig von den Eigenschaften der Matrix und den Details des iterativen L¨osungsverfahrens.
Es gibt zwei Gruppen von Verfahren
• die “klassischen” Verfahren, das Jacobi- und das Gauß-Seidel- Iterationsverfahren und darauf aufbauend Relaxationsverfah- ren, entwickelt im sp¨aten 18. Jahrhundert, werden heute noch angewandt.
• Krylov Unterraum-Methoden, die allgemeiner und oft schneller sind, jedoch ist in vielen F¨allen die Konvergenz unklar.
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) F¨ur 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ß f¨ur den Fehler.
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 f¨ur die 2-Norm oder der L¨ange des Vektors x
||x|| =
v u u t
n X 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))
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 f¨ur
||(I − A)|| < 1 oder A ist fast eine Einheitsmatrix.
(Die 2-Norm der Matrix ist die “Spektralnorm” und wird hier nicht n¨aher betrachtet, siehe Kapitel 2 zu Matrixnorm)
Jacobi Iterationsverfahren (1)
Verbesserung:
In vielen praktischen Problemen sind die Matrixeintr¨age auf der Dia- gonalen groß. Betrachte dann das modifizierte Problem
D−1Ax = D−1b mit D = diag(A) Diese ¨Anderung ¨andert die L¨osung nicht!
D−1A ≈ I bzw. ||(I− D−1A)|| < 1 ist aber eher erf¨ullt als ||(I −A)|| < 1 Beispiel:
A =
5 2 1 2 7 4
−1 2 8
, D−1 =
1/5 0 0
0 1/7 0
0 0 1/8
, I−D−1A =
0 2/7 1/8 2/5 0 4/8
−1/5 2/7 0
,
Matrizen, deren gr¨oßte Eintr¨age auf der Diagonalen sind, kommen z.B.
bei Beleuchtungsproblemen (Kapitel 2) oder bei sogenannten partiel- len Differentialgleichungen vor (Kapitel 10).
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 f¨uhrt zur selben Iterationsvorschrift:
Dx(t+1) = b + (D − A)x(t)
Jacobi Iterationsverfahren (3)
Elementweise geschrieben:
ajjx(jt+1) = bj + ajjx(jt) − X
k
ajkx(kt)
x(t+1)j = 1 ajj
bj − X
k6=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
• F¨uhre das solange durch, bis der Residuenvektor klein ist
Gauß-Seidel Verfahren (1)
Weitere Verbesserung: Zerlege die Matrix in 3 Teile:
• D: Diagonalteil von A
• −L: Unterdiagonalteil von A
• −U: Oberdiagonalteil von A L¨ose
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)
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).
Gauß-Seidel Verfahren (3)
Elementweise geschrieben:
(D − L)x(t+1) = b + Ux(t) ajjx(t+1)j + X
k<j
ajkx(t+1)k = bj − X
k>j
ajkx(t)k .
Berechnet x(1t+1), und damit x(2t+1) usw.
a11x(t+1)1 = b1 − X
k>1
a1kx(t)k a22x(t+1)2 + a21x(t+1)1 = b2 − X
k>2
a2kx(t)k a33x(3t+1) + a31x(1t+1) + a32x(2t+1) = b3 − X
k>3
a3kx(kt) ... = ...
Relaxationsverfahren
Die Wahl der Matrix, mit der die urspr¨unglich Fixpunktgleichung mo- difiziert wird, ist im Prinzip beliebig (siehe Pr¨akonditionierung)
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.
Pr¨ akonditionierung
Verallgemeinerung der obigen Methoden: Anstatt das Originalproblem zu l¨osen, betrachte
M−1Ax = M−1b
M sollte leicht zu invertieren sein und in der N¨ahe von A liegen, damit M−1A leicht zu berechnen ist und in der N¨ahe der Einheitsmatrix liegt.
Die Richardson Iterationsvorschrift f¨ur das pr¨akonditionierte 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 Pr¨akonditionierungsmatrizen M = D bzw. M = D − L.
Krylov Unterraum-Methoden (1)
Problem:
Die bis jetzt vorgestellten Methoden konvergieren nur, wenn die Itera- tionsmatrix ||C|| = ||1 − M−1A|| klein ist.
F¨ur das Jacobi-, Gauß-Seidel- und SOR-Verfahren bedeutet das, dass die Nebendiagonalelemente von A im Vergleich zu den Diagonalele- menten klein sein m¨ussen.
Krylov Unterraum-Methoden beruhen auf einer Verallgemeinerung der bisher vorgestellten Fixpunktgleichungen. Aus dem pr¨akonditio- 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 m¨ogliche Richtung annehmen kann, so dass die L¨osung gefunden werden kann, und dass der Betrag des Residuums ||b − Axk|| immer kleiner wird.
Krylov Unterraum-Methoden (2)
Was ist ein Krylov Unterraum?
Beispiel: F¨ur 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)
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 X j=t−1
r(j)
= x(0) +
t X j=0
r(j)
=
t X j=0
(I − A)j b
Krylov Unterraum-Methoden (4)
Die Iterationsl¨osung x(t) = Pt−1j=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 Iterationsl¨osung xt liegt bei einem Krylov-Unterraumverfahren in Kt(A,b).
Je nachdem, wie die Iterationsl¨osung innerhalb des Unterraums bestimmt wird, gibt es verschiedene Iterations-Strategien.
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) ∈ Kt⊥
• 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.
Einschub: Schreibweisen f¨ ur 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, h¨aufig bei Ingenieuren und in der Schule verwendet
* ha,bi 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)
CG (1)
Das konjugierte Gradientenverfahren (CG) ist Grundlage vieler moder- ner Iterationsverfahren (entwickelt 1952 von Stiefel und Hestens).
• Es konvergiert (falls keine Rundungsfehler vorliegen) f¨ur positiv definite und symmetrische Matrizen
a(x) = xTAx = (x, Ax) > 0; Ai,j = Aj,i
der Gr¨oße n × n in n Schritten und ist somit eine schnelle und sichere Methode zur L¨osung des Gleichungssystems.
• Die Rundungsfehler f¨uhren bei großen Systemen jedoch dazu, dass man in vielen F¨allen nicht n sondern ca. 3n Schritte anwendet.
• Es ist das einzige iterative Verfahren, f¨ur das die Konvergenz im Allgemeinen bewiesen wurde.
• Da jeder Schritt einer Matrixmultiplikation entspricht, lohnt sich das Verfahren besonders f¨ur d¨unn besetzte Matrizen.
CG (2)
Grundlage ist ein Optimierungsproblem. Die L¨osung von Ax = b
ist ein Minimum der Funktion f(x) = 1
2(x, Ax) − (b,x)
und umgekehrt. Begr¨undung: Sie x der L¨osungsvektor von Ax = b, so gilt f¨ur 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”.
CG (3)
Die Methode des “steepest descent” der Gradientenverfahren:
• Richtung: Beginnt man mit einem Vektor x + p und m¨ochte ent- lang des steilsten Abstiegs das Minimum erreichen, so muss man entlang 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 bestimmt f¨ur 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)
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)
• F¨uhre 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 f¨ur positiv definite und symmetri- sche Matrizen, aber meist sehr langsam, da die Richtung des steilsten Abstiegs zu einer oszillierenden Bewegung der L¨osung f¨uhren kann.
CG (5)
Beispiel aus G.Opfer:
Bestimme die L¨osung des Gleichungssystem 1 0
0 10
!
· x1 x2
!
= 0
0
!
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 0 0 10
!
· x1 x2
!
− 0 0 · x1 x2
!
= 1
2(x21 + 10x22)
CG (6)
W¨ahle als Startvektor x(0)T = (1.0, 0.1)
Der L¨osungsvektor oszilliert zur L¨osung
-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)
CG (7)
Die Methode der “konjugierten Richtung”:
Ausgangspunkt f¨ur Methode des “steepest descent” war die Funktion f. Als Richtung wurde der Gradient festgelegt und die L¨ange 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 Verst¨andnis, ohne Beweise):
• Die Richtungen, in der sich x(i) entwickelt, seien in einem Unter- raum P(i) =< p(0), p(1), . . . p(i−1) >
CG (8)
• Der Vektor x(i) l¨asst sich dann schreiben als x(i) = x(i−1) + α(i−1)p(i−1) = x(0) +
i−1 X j=0
α(j)p(j)
• W¨ahle als n¨achste Richtung zur Berechnung von x(i+1) nicht den Gradienten, sondern eine Richtung p(i), so dass der 2. Summand in f(x + αp) f¨ur alle Betr¨age von x in Richtung x(i) wegf¨allt
(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 f¨ur j ≤ i
CG (9)
• Die Vektoren p(0),p(1), . . . ,p(i) heißen paarweise zu A-konjugierte Richtungen, da gilt
(p(k), Ap(j)) = 0 f¨ur 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)
CG (10)
Algorithmus: (bis jetzt)
Verwende einen Satz von konjugierten Suchrichtungen pi. Setze x(0) und r(0) = b − Ax(0)
Dann f¨uhre 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) g¨unstig zu w¨ahlen, dass sich schnell eine gute N¨aherung ergibt.
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 gew¨ahlt, 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)) f¨ur i = 0, . . . , m−1
CG (12)
Da
(p(m−1), Ap(i)) = 0 f¨ur i < m − 1 folgt
β(m−1) = − (r(m), Ap(m−1)) (p(m−1), Ap(m−1))
Ohne Beweis: Jetzt gilt f¨ur den neuen Residuenvektor (r(i+1),p(j)) = 0 f¨ur 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 vergr¨oßert, ist sicher gestellt, dass das Verfahren nach n Schritten eine L¨osung findet.
Nach einigen Umrechnungen ergibt sich CG-Algorithmus.
CG (13)
Endg¨ultiger Algorithmus (ohne Beweis!):
Starte mit p(0) = r(0) = b−Ax(0). Dann f¨uhre folgende Schritte durch:
α(i−1) = |r(i−1)|2
(p(i−1), Ap(i−1))
x(i) = x(i−1) + α(i−1)p(i−1) r(i) = r(i−1) − α(i−1)Ap(i−1) β(i) = |r(i)|2
|r(i−1)|2
p(i) = r(i) + β(i)p(i−1)
CG (14)
Nochmal das Beispiel aus G.Opfer
1 0 0 10
!
· x1 x2
!
= 0
0
!
W¨ahle als Startvektor wieder x(0)T = (1.0,0.1) Der L¨osungsvektor oszilliert nicht mehr und die L¨osung 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)
CG (15)
Rechenaufwand pro Iteration 1 Matrixmultiplikation A,
2 Skalarprodukte, 3 AXPYs
Es werden (theoretisch) n Iterationen ben¨otigt.
Einfache Erweiterung f¨ur nicht-symmetrische Matrizen
Ist die Matrix nicht symmetrisch und positiv definite, betrachte ATAx = ATb
und konstruiere die L¨osung in Ki(ATA, ATb) unter Verdoppelung der Anzahl der Matrix-Vektor Multiplikationen (geht besser durch Erwei- terung des Krylov-Unterraums).
