6 Geconjungeerde gradienten
Laat A ∈ IRn×neen symmetrische positief definiete matrix zijn, d.w.z.
AT = A en er is een γ > 0 zodat xTA x ≥ γ xTx voor alle x ∈ IRn, (6.1) dan is het oplossen van het stelsel vergelijkingen A x = b voor gegeven b ∈ IRn equivalent met het minimaliseren van de functie
x7→ F (x) := (xb− x)TA (xb− x) , (6.2) waar xb := A−1b de oplossing is van A x = b . Omdat A (strikt) positief definiet is, is F positief als x 6=xb. De functie F kunnen we ook anders schrijven,
F (x) := (A−1b− x)TA (A−1b− x) = xTA x − 2bTx+ bTA−1b. (6.3) Omdat de derde term bTA−1b constant is heeft deze geen invloed op de argument x dat F mini- maliseert en mogen we deze term weglaten. In het vervolg gebruiken we dus de objectfunctie
F (x) := xTA x − 2bTx. (6.4)
We beschouwen de volgende iteratieve algoritme (afdalings- of descentmethode) voor het minimali- seren van F :
kies een startvector x0; r0:= b − A x0; k := 0 ; while rk6= 0 do
kies afdaalrichting pk en minimaliseer F langs de lijn xk+ λpk
d.w.z. kies λk zo dat F (xk+ λkpk) ≤ F (xk+ λpk) voor alle λ ; xk+1:= xk+ λkpk;
rk+1:= b − A xk+1 = rk− λkA pk; k := k + 1 ;
end
(6.5)
...............................................................................................................................................................................................................................
...............................................................................................................................................................................................................................................................................................................................................................................................................................................................................................
............
..................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................
...
...
...
...
...
...
...
......... ...
.................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................
...
...
....
³³
³³
³³
³³
³³
³³
³³
³³
³³
³³ )
x1
AA AA
AA AA
AAU x2
XXXXXXXX
XXXXXXXz x3
¡¡¡¡¡µ x4
Figure 12: Hoogtelijnen van de objectfunctie en opeenvolgende iteranden en zoekrichtingen in een afdaal- methode.
In deze algoritme kiezen we een beginpunt x0 en een zoekrichting p0en we zoeken het minimum van de functie F langs de lijn x0 + λp0. Langs zo’n lijn is de functie een eenvoudige parabool, waarvan we het minimum eenvoudig kunnen bepalen. In dit minimum is de functie uiteraard
kleiner dan in het startpunt. Dit gaat dan onze volgende iterand zijn, van waaruit we de procedure herhalen. In de k-de slag zoeken we dus vanuit xk het minimum van λ 7→ F (xk + λpk). Dit minimum λk is eenvoudig te bepalen door te differenti¨eren
0 = d
dλF (xk+ λpk) = 2pT(A(xk+ λpk) − b) zodat λk= pTk(b − A xk)
pTkA pk = pTkrk
pTkA pk (6.6) Als volgende iterand kiezen we dan xk+1 := xk+ λkpk. Aangezien we een oplossing van Ax = b zoeken, berekenen we ook het residu rk+1 := b − A xk+1, om te zien hoe ver we van ons doel verwijderd zijn. Voor het residu in (6.5) zijn twee (mathematisch) equivalente formules gegeven, het “echte” residu b − A xk+1 en het “recursieve” residu rk − λkA pk. In de praktijk kan het recursieve residu door afroundfouten op den duur echter sterk van het echte residu gaan afwijken.
Voor een zekere keuze van zoekrichtingen is het gedrag van een descentmethode geschetst in figuur 12. We zien, dat de lijn waarlangs we de functie minimaliseren in het minimum precies raakt aan de hoogtelijn (of niveaukromme) van de functie door dit minimum. Deze figuur suggereert, dat de methode wel zal convergeren als we de hoek tussen twee opeenvolgende zoekrichtingen niet te klein nemen. Er blijft dus een enorme keuzevrijheid over, die we kunnen gebruiken voor het optimaliseren van de methode.
Opgave 1: Bij de methode van Gauss–Seidel (relaxatiemethode) splitsen we de matrix A in de diagonaal D, het strikte linksonderstuk L en het rechtsbovenstuk U , zodat A = D + L + U . De algoritme kun je dan op de volgende manieren formuleren:
kies een startvector y0; v0:= b − U y0; k := 0 ; while vk− vk−16= 0 do
Los yk+1 op uit (D + L)yk+1= vk; vk+1= b − U yk+1; k := k + 1 ; end
kies een startvector y0; r0:= b − A y0; k := 0 ; while rk 6= 0 do
Los uk op uit (D + L)uk= rk; yk+1:= yk+ uk;
rk+1= −U uk; k := k + 1 ; end
(6.7)
Omdat D + L een onderdriehoeksmatrix is, kunnen we het stelsel (D + L)yk+1= rk eenvoudig oplossen met een voorwaartse substitutie. Waarom is D + L niet singulier?
Laat zien dat dit een descentmethode is, waarin de zoekrichtingen cyclisch de verzameling basisvectoren {ek| k = 1 · · · n} met (ej)i= δij (Kroneckers delta) doorlopen:
pk:= ej(k) met j(k) = 1 + k mod n ;
als {xk} de rij iteranden is van de descentmethode (6.5) met startvector x0:= y0, dan geldt xnj+i= (yj+1,1, · · · , yj+1,i, yj,i+1, · · · , yj,n)T
waar yj,i de i-de component is van yj.
In de gradi¨entmethode kiezen we vanuit xk als zoekrichting de richting waarin de objectfunctie F het snelst daalt. Dit is de richting van de gradi¨ent, dus
pk = − ∇F|x=xk = 2(b − A xk) = 2rk. (6.8) De algoritme heeft dan de vorm:
kies een startvector x0; r0 := b − A x0; k := 0 ; while rk6= 0 do
λk= rTkrk rTkA rk
; {minimum langs afdaalrichting rk} xk+1:= xk+ λkrk;
rk+1 := b − A xk+1= rk− λkA rk; k := k + 1 ;
end
(6.9)
De zoekrichting in k-de stap is in dit geval dus de richting van het residu rk. Vanwege de mini- malisatie in de vorige (k−1–ste) stap is de afgeleide van F in xk in de richting rk−1 = pk−1 gelijk aan nul, zodat de nieuwe zoekrichting loodrecht staat op de vorig.
Geconjungeerde zoekrichtingen. Bij de gradientmethode kiezen we, om vanuit xk de volgende benadering te berekenen, als zoekrichting p de richting van de gradient, omdat F in die richting het snelst daalt (althans in xk). De richting p is dus zo, dat
∂F
∂p|x=xk = lim
t→0
F (xk+ tp) − F (xk)
t kpk = (∇F )Tp
kpk = 2(b − Axk)Tp
kpk (6.10)
maximaal is. Hieruit volgt inderdaad, dat de maximaliserende richting gelijk is aan de gradient van F in xk en dus loodrecht staat op het niveau-oppervlak {x ∈ IRn| F (x) = F (xk)} van F door xk. Nu zijn deze niveau-oppervlakken in het algemeen ellipsoiden en geen bollen. De normaal in xkzal dus in het algemeen niet door het minimum van F gaan, tenzij xkeen heel speciaal punt is. Alleen bij een bol gaat iedere normaal door het middelpunt. Alleen als de niveauoppervlakken bollen zijn hebben we aan ´e´en iteratie genoeg om het minimum van F te bepalen. Er is ´e´en metriek waarin de niveauoppervlakken van F bollen zijn en deze is gegeven door
x7→ kxkA:=√
xTA x met bijbehorende bilineaire vorm (x , y)A7→ xTA y . (6.11) Omdat A symmetrisch en positief definiet is, is deze metriek niet gedegenereerd en is de bijbe- horende bilineaire vorm een inproduct in IRn. Uit (6.10) volgt dan dat de optimale zoekrichting t.o.v. deze metriek de richting is die de functionaal
p7→ 2(b − Axk)Tp ppTA p
maximaliseert. Deze wordt dus gegeven door p = A−1(b − Axk), maar hiervoor moeten we een stelsel van de vorm Ax = b oplossen. De steilste helling in de metriek (6.11) is dus niet te bepalen.
We weten echter wel, dat deze steilste helling ligt in het (hyper-)vlak dat in xk A-loodrecht (in de zin van 6.11) staat op de vorige zoekrichting pk−1en dat het gezochte minimum ook in dit loodvlak ligt. Dit betekent, dat we de verdere zoekactie kunnen beperken tot dit loodvlak en dus tot een deelruimte van kleinere dimensie. Als we de zoekrichting pk niet alleen A-loodrecht op de vorige zoekrichting maar op alle voorgaande zoekrichtingen kiezen, en dus als
pTkA pj = 0 , j = 0 , · · · , k − 1 , (6.12) dan wordt de dimensie van de ruimte, waarin we het minimum zoeken, in iedere slag met ´e´en ver- laagd, zodat we na n slagen het minimum gevonden hebben. Een stelsel vectoren {p0, · · · , pn−1}, dat aan relatie (6.12) voldoet, heet een A-geconjungeerd stelsel. Aangezien A symmetrisch en positief definiet is, vormt een dergelijk stelsel vectoren een basis in IRn.
Dit is het idee achter de stelling, dat algoritme (6.5) in hooguit n slagen naar de oplossing van Ax = b convergeert, als we achtereenvolgens de vectoren van het A-geconjungeerde stelsel {p0, · · · , pn−1} als zoekrichtingen kiezen. Een preciese formulering van deze stelling is de volgende:
Stelling 6.1 Zij A ∈ IRn×n een symmetrische positief definiete matrix, zij {p0, · · · , pn−1} een A-geconjungeerd stelsel vectoren en zij x0∈ IRn een willekeurige startvector voor de algoritme
r0 := b − A x0; k := 0 ; while k rkk > 0 do
λk:= pTkrk pTkA pk
; {minimum langs afdaalrichting pk} xk+1:= xk+ λkpk;
rk+1 := b − A xk+1= rk− λkA pk; k := k + 1 ;
end .
(6.13)
Dan is de laatste iterand xk de gezochte oplossing van het stelsel vergelijkingenAx = b .
Bewijs: Laat xb de oplossing zijn van Ax = b . Daar het stelsel {p0, · · · , pn−1} een basis in IRn is, kunnen wexb− x0 op unieke wijze schrijven als een lineaire combinatie van deze basisvectoren,
b
x− x0 =
n−1X
i=0
αipi en dus r0:= b − Ax0= A(xb− x0) =
n−1X
i=0
αiApi. (6.14) Bijgevolg geldt:
λ0= pT0r0
pT0A p0
=
n−1X
i=0
αipT0A pi
pT0A p0
= α0
en dus vinden we b
x− x1=
n−1X
i=1
αipi ∈ {u ∈ IRn| uTA p0= 0} en rk=
n−1X
i=k
αiA pi⊥ p0. Analoog vinden we in alle stappen van de algoritme λk= αk,
b
x− xk=
n−1X
i=k
αipi ∈ {u ∈ IRn| uTA pj = 0 , j = 0 · · · k − 1}
en rk=
n−1X
i=k
αiA pi ⊥ pj (j = 0 · · · k − 1) .
(6.15)
Het is mogelijk, dat algoritme (6.13) stopt na minder dan n stappen, maar het residu is na n stappen zeker gelijk aan nul en als het residu nul is, dan is de bijbehorende xk gelijk aan de oplossing bx.
Opmerking 6.2 : Uit (6.15) zien we datxb− xk A-loodrecht staat op alle voorgaande zoekrichtin- gen, zoals al eerder opgemerkt, en dat rk (gewoon) loodrecht staat op alle voorgaande zoekrichtin- gen. We kunnen dus ook zeggen, dat de minimalisatie van F (xk+ λpk) in de k-de slag niet alleen het minimum geeft op de lijn door xk parallel aan pk, maar zelfs in de gehele affiene deelruimte (of lineaire vari¨eteit) door x0 parallel aan vect{p0, · · · , pk} .
Met een A-geconjungeerd stelsel zoekrichtingen wordt de iteratieve methode (6.5) dus eindig.
Als we daarbij echter het gehele A-geconjungeerde stelsel vooraf zouden moeten uitrekenen en in het geheugen moeten bewaren, dan zou de methode minstens zoveel tijd en geheugenruimte nodig hebben als een direkte methode zoals Gauss-eliminatie. Uit de algoritme zien we echter dat we de zoekrichting pk pas in de k+1–ste slag nodig hebben en dat we dus de keuze kunnen uitstellen tot deze slag. De enige beperking bij deze keuze is, dat pk A-loodrecht moet staan op alle voorgaande zoekrichtingen. De vraag is nu of we zo’n zoekrichting kunnen vinden met weinig werk, en dus zonder alle orthogonalisaties expliciet uit te voeren.
In (6.15) zagen we, dat rk(gewoon) loodrecht op alle voorgaande zoekrichtingen {p0, · · · , pk−1} staat; we zullen laten zien, dat we een collectie zoekrichtingen kunnen vinden die zo is, dat rk ook A-loodrecht staat op al de voorgaande zoekrichtingen. Als we zo’n collectie hebben, dan hoeven we rk slechts A-loodrecht op pk te zetten om een A-geconjungeerde zoekrichting te vinden voor de k+1–ste slag; we kiezen dan
pk+1 := rk+1−rTk+1A pk
pTkA pk
pk (6= 0 als rk+16= 0 , omdat pk ⊥ rk+1) . (6.16) Veronderstel dat λj 6= 0 als j < k. (d.w.z. veronderstel krjk 6= 0). Uit (6.13) zien we dan
rj+1= rj− λjA pj zodat A pj = 1 λj
(rj− rj+1) . (6.17)
Voor het inproduct met rk vinden we dan rTkA pj = 1
λjrTk(rj− rj+1) . (6.18) Als rj een lineaire combinatie van {p0, · · · , pj} is, dan is het rechterlid van (6.18) nul voor j = 0 , · · · , k − 1 , omdat volgens (6.15) rTkpi = 0 voor i = 0 , · · · , k − 1 . Als we de keuze (6.16) in iedere stap hebben gedaan, is aan deze voorwaarde automatisch voldaan. Bovendien impliceert deze keuze, dat rk= 0 en xk =xb als λk = 0 , immers uit de definitie van λk in (6.13) en de keuze (6.16) zien we:
λk= pTkrk
pTkA pk = rTkrk pTkA pk
omdat het inproduct pTk−1rk = 0 volgens (6.15); λk = 0 impliceert dus dat rk nul is en dat we de oplossing hebben gevonden.
De afdalingsmethode (6.5) met zoekrichtingen gegeven door (6.16) heet de methode der geconjungeerde gradienten(Conjugate Gradients in het Engels). De algoritme luidt alsvolgt:
Kies startvector x0; r0 := b − A x0; k := 0 ; while krkk > 0 do
if k = 0 then p0 := r0 else (a)
µk:= − rTkA pk−1
pTk−1A pk−1 = k rkk2
k rk−1k2 ; pk:= rk+ µkpk−1; (b) end
λk := pTkrk
pTkA pk = k rkk2
pTkA pk; xk+1 := xk+ λkpk; (c) rk+1:= b − A xk+1 = rk− λkA pk; (d) k := k + 1 ;
end .
(6.19)
Stelling 6.3 Zij A ∈ IRn×n een symmetrische positief definiete matrix en zij x0 ∈ IRn een willekeurige startvector voor de algoritme (6.19) dan is er m ≤ n zodat rm= 0 en xm =xb. Bewijs: Een bewijs is hierboven gegeven; we zetten de elementen nog eens op een rijtje.
Als k = 0 en r0 6= 0 dan maken we in (a) de vector p0 = r0 6= 0 zodat r0 ∈ vect{p0}. In (c–d) vinden we dan λ0 = rT0r0/rT0A r0 6= 0 ,xb− x1 ⊥A{p0} en r1 ⊥ {p0} .
Als k > 0 en rk6= 0 en als (bij inductieaanname) geldt:
i. λj 6= 0 voor j = 0 · · · k − 1 ,
ii. {p0, · · · , pk−1} is een A-geconjungeerd stelsel,
iii. vect{p0, · · · , pj} = vect{r0, · · · , rj} = Kj+1(A, r0) := vect{r0, Ar0, · · · , Ajrj} voor j = 0 · · · k − 1 , waar Kk(A, r0) de k-de Krylovruimte van A en r0 genoemd wordt.
iv. xb− xk⊥A{p0, · · · , pk−1} en rk ⊥ {p0, · · · , pk−1} ,
dan volgt in (b), dat pk 6= 0 omdat rk ⊥ pk−1 en dat bij constructie pk A-loodrecht staat op pk−1. Omdat λjrTkA pj = rTk(rj+1− rj) = 0 voor j = 0 · · · k − 2 (zie 6.18), staat pk ook A-loodrecht op alle voorgaande zoekrichtingen, zodat {p0, · · · , pk} weer een A-geconjungeerd stelsel is met vect{p0, · · · , pk} = vect{r0, · · · , rk}.
In (c-d) volgt tenslotte, datxb− xk+1⊥A{p0, · · · , pk} en rk+1 ⊥ {p0, · · · , pk} . Na hoogstens n slagen is het residu nul en is de oplossing bereikt.
Opgave 2: Laat zien, dat we de getallen λk en µk in (6.19) ook kunnen berekenen met
λk := rTkrk
pTkA pk en µk:= rTkrk
rTk−1rk−1. (6.20)
Tesamen met de twee manieren om het residu (wel of niet recursief) geeft dit 8 (analytisch) equivalente manieren om de algoritme te implementeren. Bepaal voor ieder van deze manieren de hoeveelheid werk in termen van aantallen matrix-vector vermenigvuldigingen, inproducten en vector-updates (van de vorm x:= x + αy).
Deze algoritme, ge¨ıntroduceerd door Hestenes en Stiefel [1], geeft in theorie dus een eindige methode om de oplossing van Ax = b te berekenen. Deze methode is vooral geschikt voor ijle matri- ces, d.w.z. matrices waarvan de meeste elementen nul zijn zodat een matrix-vector vermenigvuldig- ing veel minder dan O(n2) flops vraagt. Helaas is de eindigheid van de algoritme niet bestand tegen de eindige precisie van een computer. Door afrondfouten staat de berekende vector pk niet exact A-loodrecht op al zijn voorgangers. De afwijking t.o.v. de loodrechte stand tussen pk en pj wordt groter naarmate het verschil | k − j | groter wordt.
Figure 13: Het trampolinerooster met 8 knopen horizontaal en 6 verticaal. De randknopen zijn vast. Rond knoop (4,3) is het gebied geschetst waarvan de totale massa op deze knoop drukt.
Voorbeeld: We willen de vorm van een trampoline (of bedspiraal) met afmetingen ℓ×b berekenen, als we deze belasten met een gewicht g(x, y) (per oppervlakte-eenheid). Modelleer de trampoline als een rechthoekig array van m × n knopen, verbonden door veren van lengte h, zie fig. 13, zodat dus ℓ = mh en b = nh . We kunnen de verticale kracht Fi,j op knoop (i, j) schrijven als de som van de verticale krachten langs de vier veren. Deze krachten zijn evenredig met het hoogteverschil, zodat
Fi,j = S(ui,j−1− ui,j) + S(ui,j+1− ui,j) + S(ui−1,j− ui,j+ S(ui+1,j − ui,j) ,
als ui,j de verticale uitwijking is in knoop (i, j), S de veerconstante is en als de verschillen in de uitwijkingen klein zijn t.o.v. h. Omdat de verticale kracht Fi,j op knoop (i, j) evenredig is met het gewicht dat op een elementair vierkantje drukt en dus (ongeveer) evenredig is met h2 maal g(ih, jh) vinden we de (benaderende) vergelijking
ui,j−1+ ui,j+1+ ui−1,j+ ui+1,j− 4ui,j = h2gi,j
S met 0 < i < m en 0 < j < n . (6.21) De rand van de trampoline zit vast, zodat
u0,j = um,j = ui,0= ui,n= 0.
In de andere punten van het rooster vinden we (n − 1)(m − 1) vergelijkingen voor evenveel onbek- enden. We ordenen de (niettriviale) onbekenden en de bijbehorende rechterleden in vectoren van lengte (m − 1)(n − 1), en stellen de bijbehorende matrix op. Ga na, dat in het geval (m, n) = (5, 4) en “lexicografische” ordening van de knopen (begin linksonder en doorloop eerst alle knopen met dezelfde y-waarde) de matrix de volgende vorm heeft:
∗ ∗ ∗
∗ ∗ ∗ ∗
∗ ∗ ∗ ∗
∗ ∗ ∗
∗ ∗ ∗ ∗
∗ ∗ ∗ ∗ ∗
∗ ∗ ∗ ∗ ∗
∗ ∗ ∗ ∗
∗ ∗ ∗
∗ ∗ ∗ ∗
∗ ∗ ∗ ∗
∗ ∗ ∗
Dit is een typisch voorbeeld van een ijle matrix. Per rij zijn er hoogstens vijf elementen ongelijk aan nul, zodat een matrix-vectorvermenigvuldiging hoogstens 5(m − 1)(n − 1) flops vraagt als we voor deze matrix-vector vermenigvuldiging een routine schrijven die rekening houdt met de speciale vorm.
Dit voorbeeld kan worden opgelost met de methode de geconjungeerde gradienten. In figuur 14 zijn de residunormen getekend als functie van de iteratie-index. We zien dat het residu al tot de machineprecisie is gereduceerd lang voor het theoretische einde van het proces. We zien ook dat het echte residu (zoals verwacht) rond de machineprecisie blijft hangen terwijl het recursieve residu gewoon verder daalt en kennelijk geen relatie meer heeft met het echte residu. Ook zien we dat de A-orthogonaliteit van p0 en pk met het klimmen van k volledig verdwijnt.
De conclusie die we hieruit kunnen trekken is, dat CG niet moet worden gebruikt als direkte methode, maar als iteratieve, die na een aantal slagen, dat veel kleiner is dan de dimensie van het probleem, al een goede benadering van de oplossing geeft. Het was Reid [3] die als eerste in 1971 hierop wees.
Opgave 3: Laat Uk het k-de Chebyshev polynoom van tweede soort zijn (zie syllabus 7.c opgave 3). Laat zien dat de vector met componenten uk,j := Uk−1(ξ) Uj−1(η) een eigenvector is van de matrix in het linkerlid van (6.21) behorend bij de eigenwaarde 2ξ + 2η − 4 , als ξ een nulpunt is van Um−1 en η een nulpunt van Un−1, en dat bijgevolg het conditiegetal (t.o.v. de Euklidische norm) van de matrix gelijk is aan
κ2= 2 + cosmπ + cosπn 2 − cosmπ − cosπn
≈ 4 n2
π2 als n = m (6.22)
Geconjungeerde gradienten als iteratieve methode. Om Geconjungeerde gradienten te kun- nen vergelijken met andere iteratieve methoden herschrijven we (6.19) door de vectoren pk te elimineren met gebruik van (6.17):
rk+1= rk− λkA pk= rk− λkA(rk+ µkpk−1) = rk− λkA rk+λkµk
λk−1(rk− rk−1) zodat we effectief de drietermsrecursierelatie vinden:
rk+1= µ
1 +λkµk
λk−1 − λkA
¶
rk− λkµk λk−1
rk−1. (6.23)
Als we de rij polynomen {pk} defini¨eren door de drietermsrecursie pk+1(x) :=
µ
1 −λkµk λk−1 − λkx
¶
pk(x) + λkµk
λk−1 pk−1(x) (k ≥ 2), p0(x) := 1 en p1(x) := 1 − λ0x ,
(6.24)
0 20 40 60 80 100 120 140 160 180 200 10-20
10-15 10-10 10-5 100
residu-norm van Jacobi en CG voor nxn vierkant, n = 41
iteratie index
residu norm
Jacobi
echt cos v/d hoek tussen p(k) en p(0) bij echt residu
cos v/d hoek tussen p(k) en p(0) bij recursief residu
Chebyshev bovengrens
Figure 14: Het oplosproces voor het stelsel vergelijkingen (6.21) met n = m = 41 zodat de dimensie van de oplosruimte 1600 is. Als functie van de iteratie-index zijn uitgezet: de norm van het residu van Jacobi- iteratie en de normen van de residuen van geconjungeerde gradienten met echt met recursief residu. Voor beide varianten is ook de absolute waarde van de cosinus van de A-hoek tussen p0 en pk uitgezet.
dan geldt voor iedere k (ga na!):
rk= pk(A) r0 en pk(0) = 1 . (6.25) Op dezelfde manier elimineren we pk uit xk:
xk+1 = xk+ λkpk= xk+ λk(rk+ µkpk−1) = xk+ λk(b − A xk) +λkµk
λk−1(xk− xk−1) , xk+1−bx = xk−xb− λk(A xk− Ax) +b λkµk
λk−1((xk−xb) − (xk−1−bx)) = pk(A)(x0−x)b
= x0−xb+ (pk(A) − 1)A−1A(x0−x) = xb 0−xb+ (1 − pk(A))A−1r0. We vinden zo een polynoom qk van graad k−1 waarvoor geldt:
xk+1 = x0+ qk(A) r0 met qk(x) := 1 − pk(x)
x . (6.26)
We herschrijven de CG-algoritme (6.19) hiermee formeel alsvolgt:
kies een startvector x0; r0 := b − A x0; k := 0 ; while rk6= 0 do
xk+1 := x0+ qk(A)r0; rk+1 := pk+1(A)r0; k := k + 1 ;
end
(6.27)
We zien uit (6.25) dat rk en xk+1− x0 elementen zijn van de k-de Krylov-ruimte Sk van A en r0
die opgespannen wordt door de vectoren r0 · · · Akr0,
Sk:= vect{r0, A r0, A2r0, · · · , Akr0} (6.28) Ter vergelijking beschouwen we een ander proces in diezelfde Krylov-ruimte, successieve sub- stitutie. We schrijven A x = b als x = x + b − A x , we kiezen een startvector x0 en de iteratie xn+1 = xn+ b − A xn. Het residu is dan rn := b − A xn en rn+1− rn = A xn+1− A xn = A rn. Alles bijeen vinden we dus:
kies een startvector x0; r0 := b − A x0; k := 0 ; while rk6= 0 do
xk+1 := xk+ rk; rk+1 := (1 − A)rk; k := k + 1 ;
end
(6.29)
We zien, dat rk:= (1 − A)kr0= pk(A)r0 en xk = x0+ r0+ · · · + rk−1 = x0+
k−1X
j=0
(1 − A)jr0 = x0+ (1 − (1 − A)k)A−1r0= x0+ qk(A)r0, met pk(x) := (1 − x)k en qk(x) := (1 − pk(x))/x . Deze methode heeft zo dus dezelfde vorm als (6.27). Convergentie treedt op als limk→∞k(1 − A)kk = 0 in een of andere matrixnorm, d.w.z. als de absolute waarden van alle eigenwaarden van A strikt kleiner dan 1 zijn.
Analoog aan Cesaro-sommatie kunnen we de door (6.29) voortgebrachte rij {xk} omzetten in een nieuwe rij {yk} , die sneller naarxb convergeert, door voor yk een geschikte lineaire combinatie van {x0 · · · xk} te nemen,
yk:=
Xk j=0
γkjxj met
Xk j=0
γkj = 1 .
Als x0 = x, dan geldt xb k = xb voor alle k; wegens de voorwaarde Pkj=0γkj = 1 geldt dan ook yk=xb. Voor het residu sk:= b − A yk betekent de voorwaarde op de som:
sk:= b − Xk j=0
γkjA xj = Xk j=0
γkjA(bx− xj) = Xk j=0
γkjrj = Xk j=0
γkj(1 − A)jr0 =: πk(A)r0, waar πk:=Pkj=1γkjpj opnieuw een polynoom van graad k is dat voldoet aan πk(0) = 1 voor alle k vanwege de somconditie Pkj=0γkj = 1 . We vinden analoog aan (6.26) het geassocieerde polynoom ϕk waarvoor yk voldoet aan de relatie
yk= x0+ Xk j=0
γkj(xj− x0) = x0+ Xk j=1
γkjqj(A)r0 = x0+ ϕk(A)r0 als ϕk(x) := 1 − πk(x)
x .
Zonder referentie naar de oorspronkelijke rij {xk} kunnen we de recursie voor de nieuwe rij {yk} dan analoog aan (6.27) herschrijven als:
kies een startvector y0; r0 := b − A y0; s0:= r0; k := 0 ; while rk6= 0 do
yk+1 := y0+ ϕk(A)r0; sk+1:= πk+1(A)r0; k := k + 1 ;
end
(6.30)
De enige eis voor convergentie van (6.30) is: limk→∞kπk(A)k = 0 .
Voor iedere matrix A is er zo’n rij polynomen te vinden, kies bijvoorbeeld alle pk met k ≥ n gelijk aan het karakteristieke polynoom ΠAvan A, dan geldt automatisch kpk(A)k = 0 (k ≥ n), zie (6.31).
De constructie van het karakteristieke polynoom van A (voor grote n) en in het algemeen ook van een rij polynomen waarvoor (6.30) convergent is, is echter geen eenvoudige opgave.
Voor het geval dat A symmetrisch en positief definiet is, hebben we echter zo’n methode gevon- den, nl. geconjungeerde gradienten (6.19). Zoals in (6.27) aangetoond, construeert deze methode (impliciet) een rij polynomen {pk}, waarvoor het schema (6.30) convergeert in eindig veel stappen.
Bovendien is deze methode optimaal. De geconjungeerde-gradientenmethode kiest in de k-de slag het polynoom pkzo, dat de functionaal x 7→ F (x) geminimaliseerd wordt in de ruimte x0+Kk(A, r0) (zie opmerking 6.2). Omdat voor het residu van CG geldt, dat rk = pk(A)r0 ∈ Kk+1(A, r0) en rk ⊥ Kk(A, r0) (zie 6.15), is dit equivalent met minimalisatie van ̺ 7→ k̺(A)r0k over alle poly- nomen ̺ van graad k + 1 met ̺(0) = 1 . Het residu van CG is in iedere stap dus kleiner dan het residu verkregen met een andere iteratieve methoden van de vorm (6.30).
Vegelijking CG met Chebyshev iteratie: We willen de convergentiesnelheid van CG schatten, dus we wensen een (goede) bovengrens te vinden voor de norm van het residu rk = pk(A)r0 van CG.
Als A een symmetrische matrix is dan heeft deze een eigenwaardeontbinding
A = U Λ U−1, U = (u1| · · · | un) en A uk= λkuk (6.31) waarin Λ = diag(λ1, · · · , λn) een diagonaalmatrix is bestaande uit de eigenwaarden van A en waarin U een orthogonale matrix is, waarvan de kolommen de eigenvectoren zijn. Als p een poly- noom is, dan geldt
p(A) = U p(Λ) U−1 met p(Λ) = diag(p(λ1) , · · · , p(λn)) .
Als r0 = Pni=1αiui, dan is kp(A)r0k2 = Pni=1α2ip(λi)2. Een bovengrens voor deze norm hangt dus uitsluitend af van de de waarden van het polynoom p op de eigenwaarden (het spectrum) van A. Om een bovengrens voor de norm van het residu pk(A)r0 in de CG-methode te vinden zouden we de eigenwaarden van A moeten kennen en de waarden van pk op deze eigenwaarden. Dat is onbegonnen werk.
We weten echter wel, dat CG optimaal is en dat dus iedere andere serie polynomen een grotere bovengrens geeft. Bovendien weten we, dat A positief definiet is, zodat haar eigenwaarden in een interval 0 < a ≤ λj(A) ≤ b liggen. We kunnen dan een bovengrens vinden met een rij polynomen die uniform klein zijn op het interval [a, b]. De optimale polynomen hiervoor zijn de Chebyshev polynomen. Het k-de Chebyshev polynoom is gedefinieerd door
Tk(cos t) := cos kt (6.32)
zodat Tk voldoet aan de recurrente betrekking
T0 = 1 , T1(x) = x en Tk+1(x) + Tk−1(x) = 2xTk(x) voor k > 0 . (6.33) De functie t 7→ cos kt heeft op het interval [0, π] precies k + 1 maxima en minima met waarden om en om +1 en −1. Aangezien de afbeelding t 7→ cos t het interval [0, π] een-eenduidig op [−1, 1]
afbeeldt, heeft Tk dezelfde eigenschap op [−1, 1].
Voor Tk geldt de minimax eigenschap:
Stelling 6.4 Als P een polynoom van graad≤ k is met P (µ) = Tk(µ) voor een zekere µ , | µ | > 1 , dan geldt
−1≤x≤1max | P (x) | ≥ max
−1≤x≤1 | Tk(x) | ≥ 1 . (6.34)
Bewijs: Stel | P (x) | ≤ γ < 1 , dan kruist de grafiek van P die van Tk minstens k maal, omdat Tk k maal van +1 naar −1 gaat en terug en P tussen −γ en γ blijft. P − Tk heeft dus (minstens) k nulpunten binnen het open interval (−1, 1) en ook nog een er buiten (nl. µ). Omdat het een polynoom van graad k is, moet deze identiek nul zijn.
Het optimale polynoom voor een uniform kleine bovengrens op [a, b] is dus
̺k(x) := Tk(a + b − 2x
b − a )/Tk(a + b
b − a) (6.35)
Voor een bovengrens op [a, b] bewijzen we het volgende lemma:
Lemma 6.5 Als x > 1 , dan geldt:
1
2(x +px2− 1)k≤ Tk(x) ≤ (x +px2− 1)k. (6.36) Bewijs: De oplossing van de recurrente betrekking (6.33) heeft de vorm Tk(x) = αλk1(x) + βλk2(x) waar λ1,2 de wortels zijn van de karakteristieke vergelijking λ2− 2xλ + 1 = 0 , zodat
λ1 = x +px2− 1 en λ2 = x −px2− 1 . Omdat
1 = T0= α + β en x = T1(x) = α(x +px2− 1) + β(x −px2− 1) volgt α = β = 12. Bijgevolg vinden we
1
2(x +px2− 1)k≤ Tk(x) = 12(x +px2− 1)k+12(x −px2− 1)k≤ (x +px2− 1)k. De teller in ̺k is begrensd door 1 als x ∈ [a, b]. Met behulp van lemma 6.5 vinden we dus de schatting
a≤x≤bmax | ̺k(x) | ≤ 2
b + a b − a+
sµb + a b − a
¶2
− 1
−k
= 2
Ã1 −pa/b 1 +pa/b
!k
(6.37) Als we a = λmin en b = λmax kiezen dan is κ := b/a het conditiegetal van de matrix (t.o.v. de Euclidische norm). Zo vinden we tenslotte:
Stelling 6.6 Voor het residu van de geconjungeerde-gradientenmethode (6.19), geldt de volgende schatting:
krkk ≤ 2 kr0k
Ã1 −p1/κ 1 +p1/κ
!k
. (6.38)
waarκ het conditiegetal van A is.
De Ritz-waarden We keren nu terug naar de drieterms-recursierelatie (6.23) en formuleren deze als
A r0 = β0r0+ α1r1 en A rk= αk+1rk+1+ βkrk+ γkrk−1 (k ≥ 1) (6.39) met
β0 = −α1 := 1 λ0
, αk+1:= − 1
λk, γk := − µk
λk−1 en βk:= −αk+1− γk, voor k ≥ 1.
Als we alle residuen in de matrix Rk := (r0| r1| · · · | rk−1) plaatsen en met Tk de volgende tridia- gonale matrix aanduiden,
Tk:=
β0 γ1
α1 β1 γ2
α2 . .. ...
. .. ... . ..
. .. . .. γk−1
αk−1 βk−1
, (6.40)
vinden we de relaties
A Rk = RkTk+ αkrkeTk+1 en RTk A Rk= Tk. (6.41) De tweede relatie volgt uit het feit, dat rk loodrecht op de kolommen van de matrix Rk staat.
Tk is dus de restrictie van A tot de tot de k-de Krylov-ruimte opgespannen door de residuen {r0, · · · , rk−1} , te noteren met Kk(A; r0) . Deze residuen vormen tesamen een orthogonale basis.
De eigenwaarden van Tk heten de Ritz-waarden van A met betrekking tot deze Krylov-ruimte.
Opgave 4: Ga na dat de hoofddiagonaal van Tk positief is, dat de nevendiagonalen negatief zijn en dat Tk dus altijd omgezet kan worden in een symmetrische tridiagonale matrix door vermenigvuldiging met een diagonale matrix van links en met de inverse ervan van rechts.
Laat ook zien dat voor alle k ≥ 1 de volgende insluiting geldt van minimale en maximale eigenwaarden:
λmin(A) ≤ λmin(Tk) ≤ λmin(Tk−1) ≤ λmax(Tk−1) ≤ λmax(Tk) ≤ λmax(A) . (6.42)
Aangezien ook {r0, · · · , Ak−1r0} ook een basis vormt van de k-de Krylovruimte (ga na!), be- vat deze ruimte voor grote k dus een goede benadering van de eigenvector(en) van de grootste eigenwaarde(n) van A (en ook van de kleinste!). Als k voldoend groot is (maar nog veel kleiner dan n), zullen de grootste en kleinste eigenwaarden van Tk dus goede benaderingen geven van de grootste en kleinste eigenwaarden van A. De eigenwaarden van Tk kunnen snel bepaald worden met QR-iteratie, Jacobi-iteratie of met speciale met bisectiemethoden voor (symmetrische) tridi- agonale matrices. Deze combinatie van tridiagonalisatie in een Krylovruimte en de bepaling van de eigenwaarden van de (benaderende) tridiagonaalmatrix heet de algoritme van Lanczos [2]. Het zal duidelijk zijn dat ook deze methode in de praktijk zwaar te lijden heeft onder de opbouw van afrond- fouten. Desondanks kan de methode betrouwbaar ge¨ımplementeerd worden voor de benadering van de extreme (grote en kleine) eigenwaarden van een symmetrische matrix.
Als de eigenvectoren van de grootste en kleinste eigenwaarden van A goed benaderd worden in de Krylovruimte, mag je verwachten, dat de componenten van de CG-residu’s in deze richtingen ook klein zullen zijn. Dus ook het bijbehorende polynoom pk in (6.24) zal klein zijn op deze delen van het spectrum van A. Dat betekent, dat de convergentie op den duur veel sneller zal gaan dan de (pessimistische) bovengrens (6.38) verkregen door vergelijking met Chebyshev iteratie. We zien dit ook in fig. 14, waar de rechte met de kleinste helling de theoretische bovengrens voor de snelheid van Chebyshev iteratie geeft (met κ2= λmax/λmin) en waarbij de andere rechte de helling geeft op grond van eenκ die het quotient is van de op een na grootste en de op een na kleinste eigenwaardene van A . We zien in de figuur, dat het residu inderdaad steeds sneller afneemt, naarmate k groter wordt.
Preconditionering We kunnen proberen de CG-iteratie te versnellen door de (uiterste) eigen- waarden van A dichter bij elkaar (en verder van nul) te brengen door A met een geschikte matrix te vermenigvuldigen. We kunnen een matrix P zoeken, een preconditioner genaamd, zodat
λmax(P AP )/λmin(P AP ) ≪ λmax(A)/λmin(A) .
Het resultaat is dan weer een symmetrische matrix, waarop we de CG-algoritme kunnen toepassen als tevoren. Eenvoudiger wordt het echter, als we bedenken, dat we tot nu toe het standaard inproduct xT ygebruikt hebben zonder enige specifieke eigenschap ervan te gebruiken. We hadden evengoed het inproduct hx , yi := xT K y met een geschikte symmetrische positief definiete matrix K. Het is eenvoudig te verifi¨eren, dat K−1A symmetrisch en positief definiet is t.o.v. dit nieuwe inproduct. We kunnen de reeds afgeleide CG-algoritme dus geheel volgen voor het oplossen van het gepreconditioneerde systeem K−1Ax = K−1b als we het nieuwe inproduct gebruiken.
Preconditionering is voor het eerst beschreven in [4]. M&VdV kozen hiervoor een zogenaamde
“incomplete Cholesky ontbinding” van A. Hierbij worden een onderdriehoeksmatrix L (met Lkk= 1) en een diagonaalmatrix D gemaakt zo, dat L+LT hetzelfde ijlheidspatroon heeft als A, dwz. Lij 6= 0 iff Aij 6= 0 , en zo, dat (L D LT)ij = Aij voor alle (i, j) waarvoor Aij 6= 0 . Het idee hierachter is, dat deze incomplete Cholesky-factoren gelijken op de echte factoren en dat we zo een gemakkelijk te berekenen benadering van de inverse van A verkrijgen (A−1is tenslotte de beste preconditioner).
References
[1] M. Hestenes & E. Stiefel, Methods of conjugate gradients for solving linear systems, J. Research NBS, 49, pp. 409 – 436, 1952.
[2] C. Lanczos, An iteration method for the solution of the eigenvalue problem of linear differential and integral operators, J. Research NBS, 45, pp. 255 – 282, 1950.
[3] J.K. Reid, On the method of conjugate gradients for the solution of large sparse systems of linear equa- tions, Proc. Conf. on Large Sparse Sets of Linear Equations, Academic Press, New York, 1971.
[4] J.A. Meijerink and H.A. van der Vorst, An iterative solution method for linear systems of which the coefficient matrix is a symmetric M-matrix, Math.of Comp., 31, pp. 148 – 162, 1977.
[5] G.H. Golub & C.F. Van Loan, Matrix Computations, The Johns Hopkins University Press, Baltimore, Maryland, USA, 1ste druk, 1983, 2dedruk, 1988, 3de druk, 1995.
[6] R. Bulirsch & J. Stoer, Introduction to Numerical Analysis, Springer Verlag, Berlin, 1977. (Ook verkri- jgbaar in een goedkope duitstalige pocketeditie).
[7] D. Kincaid & W. Cheney, Numerical Analysis, Brooks & Cole Publishing Company, Pacific Grove, California, USA, 1991; 2de druk, 1996.
