7 Stelsels niet-lineaire vergelijkingen en minimalisatieproblemen
7.a Probleemstelling in ´e´en dimensie
Bepaal de oplossing van de volgende twee problemen:
f (x) = 0 f : (a, b) → IR voldoend glad
x∈(a,b)min F (x) F : (a, b) → IR voldoend glad (7.1) aangenomen dat zo’n oplossing bestaat. Het tweede probleem herleidt zich tot het eerste door F′(x) = 0 op te lossen, maar numeriek behoeft dit niet zinvol te zijn als de afgeleide zich niet of moeilijk laat berekenen.
7.b Intervalhalvering of binaire search
Als de tekens van f (a) en f (b) verschillen, ligt er (minstens) ´e´en nulpunt tussen en kunnen we dit vinden met de algoritme
if f(a) ∗ f(b) < 0 then repeat c:= (b + a)/2 ;
if teken(f(c)) = teken (f(a)) then a := c else b := c end until |b − a| voldoend klein .
Convergentie is verzekerd, maar vrij traag. Geen generalisatie mogelijk naar verscheidene dimensies.
7.c Successieve substitutie
Herformuleer probleem (7.1) als volgt:
x = ϕ(x) := x − f(x) (7.2)
dan kunnen we het iteratieve proces bekijken:
kies x0;
for k := 1 to · · · do xk:= ϕ(xk−1) (7.3)
Stelling 1: Als α een oplossing is van de vergelijking α = ϕ(α) en als
|ϕ′(α)| < 1 (ϕ′ continu) (7.4)
dan is er een d > 0 zodat het proces (7.3) convergeert naar α voor alle x0∈ [α − d, α + d].
Bewijs: Op grond van (7.4) is er een γ < 1 en een d > 0 zo, dat
|ϕ′(x)| ≤ γ ∀ x ∈ [α − d, α + d]
Bijgevolg geldt: als |xk− α| ≤ d, dan ook
|xk+1− α| = |f(xk) − f(α)| = |f′(ξ)| |xk− α| ≤ γ|xk− α| ≤ γd zodat |xk− α| ≤ γkd → 0 als k → ∞.
Convergentie is volgens deze stelling alleen verzekerd in een (kleine) omgeving van α als de absolute waarde van de afgeleide kleiner dan 1 is. De convergentie kan heel traag zijn (als de afgeleide dicht bij +1 of −1 is) maar ook heel snel (als de afgeleide ongeveer nul is). Het proces is onmiddellijk te generaliseren naar verscheidene dimensies:
Stelling 2: Als ϕ : D ⊂ IRn → IRn een continu differentieerbare funktie op een open gebied D ⊂ IRn is met a = ϕ(a) voor zekere a ∈ D en als
k (∇ϕ)(a) k < 1 (∇ϕ is de functionaalmatrix (∂ϕi
∂xj
)ni,j=1) (7.5) dan is er een d > 0 zo, dat B := {x ∈ IRn| k x − a k≤ d} ⊂ D en zodat het proces xn+1= ϕ(xn) convergeert naar a voor alle x0 ∈ B.
Bewijs: zie boven; vervang absolute waarden door normen.
In het successieve-substitutieproces (7.3) zit een grote vrijheid ingebakken. Omdat f (α) = 0 kunnen we i.p.v. (7.2) ook kiezen
ϕ(x) = x − f (x)
g(x) (7.6)
met een willekeurige funktie g waarvoor g(α) 6= 0. Kies bijvoorbeeld g(x) = f′(α) (mits deze constante 6= 0), dan geldt
ϕ′(x) = 1 − f′(x)
f′(α) zodat ϕ′(α) = 0
We krijgen zo een zeer snelle convergentie maar het grote probleem is natuurlijk, dat α en dus ook f′(α) onbekend zijn. Ga na, dat g(x) = f′(x) wel een goede praktische keuze is.
-0.5 0 0.5 1 1.5
-1 -0.5 0 0.5 1
f(x)
f(0) + x f’(0)
x0=0 x1
Figure 15: De funktie f met zijn raaklijn in x = 0. Deze raaklijn is een lineaire benadering van f in een omgeving van x = 0. Het nulpunt ervan geeft een betere benadering van het nulpunt van f .
7.d Newton-Raphson
Als de funktie f , waarvan we het nulpunt α willen bepalen, een continue tweede afgeleide heeft in een relevant gebied rond het nulpunt, dan kunnen we f in en (naburig punt) x0 ontwikkelen,
f (x) = f (x0) + (x − x0)f′(x0) +1
2f′′(ξ)(x − x0)2.
Voor de oplossing α van f (α) = 0 geldt dus α = x0− f (x0)
f′(x0) − 1 2
f′′(ξ)
f′(x0) (α − x0)2. (7.7) Als x0 al een goede benadering was van α, dan zal x0− f(x0)/f′(x0) een betere zijn, immers dit punt is het nulpunt van de lineaire benadering f (x0) + (x − x0)f′(x0) van f (x) in de buurt van x0; zie ook fig. 15. Zo vinden we het Newton-Raphson proces
kies x0;
for k := 0 to · · · do xk+1 := xk− f (xk)
f′(xk) end. (7.8)
Stelling 3: Als f tweemaal continu differentieerbaar is in een omgeving van α en als f′(α) 6= 0, dan convergeert het Newton-Raphson proces (7.8) lokaal kwadratisch.
Bewijs: Zij d zo, dat
maxx∈[α−d,α+d] | f′′(x) |
minx∈[α−d,α+d] | f′(x) | ≤ 1 d dan volgt uit (7.7): als xk ∈ [α − d, α + d], dan geldt
| xk+1− α |= 1 2
¯
¯
¯
¯ f′′(ξk) f′(xk)
¯
¯
¯
¯
(xk− α)2 < 1
2|xk− α |, (7.9)
zodat xk+1 ∈ [α − d, α + d] en de rij convergeert.
Opmerking: Een proces xk:= ϕ(xk−1) met ϕ(α) = α heet lokaal convergent van orde p ≥ 1 in een omgeving van α, als voor alle startpunten x0voldoend dicht bij α geldt: |xk+1−α| ≤ C |xk−α|p ∀k voor zekere C > 0. Als p > 1, dan is het proces altijd lokaal convergent (ga na!).
Een stelling zoals bovenstaande is weinig praktisch omdat het nulpunt α onbekend is. Een meer praktisch resultaat vinden we uit het volgende idee. Als f in een punt x0 voldoend klein is, als de helling van f niet al te klein is en als de kromming van f (d.w.z. de tweede afgeleide) niet al te groot is, dan kan f een doorgang door nul in de buurt van x0 niet vermijden en is het Newtonproces, dat start in x0, convergent. De preciese formulering is als volgt:
Stelling 4: Laat I ⊂ IR een open interval zijn, laat f een differentieerbare funktie zijn op I met een Lipschitz continue afgeleide en laat x0 ∈ I een punt zijn, waarvoor geldt:
(a) | f′(x) − f′(y) | ≤ γ | x − y | ∀x, y ∈ I,
(b) | 1/f′(x) | ≤ β ∀x ∈ I,
(c) | f(x0)/f′(x0) | ≤ α .
(7.10)
Als d := 1
2 α β γ < 1 en als r := α /(1 − d) zo is dat (x0− r, x0+ r) ⊂ I, dan convergeert het Newtonproces (7.8) met beginpunt x0 naar een (unieke) limiet a ∈ [x0− r, x0+ r], dat een nulpunt is van f .
Bewijs: Vanwege (7.10c) geldt:
| x1− x0 |=| f(x0)/f′(x0) | ≤ α < r (7.11) zodat x1 ∈ (x0− r, x0+ r). Vanwege (7.10a,b) geldt
|xk+1− xk| ≤ |f(xk)/f′(xk)| ≤ β |f(xk)| = β |f(xk) − f(xk−1) − f′(xk−1)(xk− xk−1)|
(7.12)
Om deze expressie te schatten gebruiken we de identiteit f (y) = f (x) +
Z y x
f′(t)dt = f (x) + (y − x)f′(x) + Z y
x
(f′(t) − f′(x))dt zodat uit (7.10a) voor alle x, y ∈ I volgt
|f(y) − f(x) − (y − x)f′(x)| ≤ γ Z y
x |t − x| dt ≤ 1
2 γ(x − y)2 Vullen we dit in formule (7.12) in, dan vinden we
|xk+1− xk| ≤ 1
2 βγ(xk− xk−1)2 (7.13)
onder de voorwaarde dat xk en xk−1 in I liggen. Dit laatste bewijzen we met induktie. Als alle iteranden {x1, x2, . . . , xk} in I liggen, dan volgt uit (7.11) en (7.13):
|x2− x1| ≤ 1
2 β γ |x1− x0|2 ≤ 1
2 βγα2= αd
|x3− x2| ≤ 1
2 β γ |x2− x1| ≤ 1
2 β γ (1
2 α2 β γ)2 = α (1
2 α β γ)3= αd3≤ αd2 en in het algemeen
|xk+1− xk| ≤ 1
2 β γ |xk− xk−1| ≤ α (1
2 α β γ)2k−1 ≤ α dk zodat
|xk− x0| ≤ |x1− x0| + |x2− xk| + · · · + |xk− xk−1| ≤ α (1 + d + d2+ · · · + dk−1) ≤ r . Bijgevolg liggen alle iteranden in het (open) interval (x0− r, x0+ r) en is de rij een Cauchy-rij met limiet in het (gesloten) interval [x0− r, x0+ r].
Een sterkere versie van deze stelling staat bekend als de stelling van Newton-Kantorowitz.
Hierbij is de begrenzing (7.10b) op de afgeleide slechts vereist in het startpunt x0 en het conver- gentiegebied wat groter. Een preciese formulering zullen we geven in de volgende paragraaf.
7.e Problemen in verscheidene dimensies
Zoals gezegd, is generalisatie van interhalvering in verscheidene dimensies niet mogelijk, terwijl de generalisatie van successieve substitutie en van Newton(achtige) processen heel gemakkelijk is. We zijn vooral ge¨ınteresseerd in Newton(achtige) processen.
Laten er n funkties f1· · · fngegevens zijn van n onbekenden x1· · · xn, waarvoor we een gemeen- schappelijk nulpunt willen vinden. We noteren de funkties en de onbekenden als vektoren in IRn. Gegeven is dus een (open) gebied D ⊂ IRn (b.v. een bol) waarop de funktie f : D → IR gedefinieerd is. We nemen aan dat alle componenten van f continue (parti¨ele) tweede afgeleiden hebben. Definieer ∇f (de Jacobiaan) als de matrix van parti¨ele afgeleiden van f,
(∇f)ij := ∂fi
∂xj
, (7.14)
dan is er (volgens de Taylorontwikkeling) een konstante γ, zo dat
kf(x) − f(y) − ∇f(y)(x − y)k ≤ γkx − yk2. (7.15)
Analoog aan het eendimensionale geval zal het Newton-Raphson proces,
xn+1= xn− (∇f(xn))−1f(xn) , (7.16) lokaal convergent zijn in een omgeving van een nulpunt a als de Jacobiaan ∇f in dit punt inver- teerbaar is:
Stelling 5: Als f een nulpunt a heeft, als de Jacobiaan ∇f(a) in dat punt inverteerbaar is en als alle parti¨ele tweede afgeleiden continu zijn in een omgeving van a, dan is er een d > 0 zo dat het Newtonproces (7.16) convergeert naar a voor iedere startpunt x0 met kx0− ak ≤ d .
Bewijs: Volledig analoog met het eendimensionale geval.
Op dezelfde manier kan stelling 4 vertaald worden naar verscheidene dimensies en bewezen worden. Een iets sterkere variant ervan is de volgende stelling van Kantorowitz. Hierin gebruiken we voor de (open) bol met straal d om y ∈ IRn de volgende notatie:
S(y, d) := {x ∈ IRn| kx − yk < d} .
Stelling 6: (Newton-Kantorowitz) Laat f een continu differentieerbare afbeelding zijn van een open gebied D ⊂ IRn naar IRn die voldoet aan de volgende eigenschappen
(a) k∇f(x) − ∇f(y)k ≤ γkx − yk ∀x, y ∈ D
(b) k(∇f(x0))−1f (x0)k ≤ α voor zekere x0∈ D (c) k(∇f(x0))−1k ≤ β
(7.17)
met α, β en γ zo, dat h := αβγ ≤ 1
2 en S(x0, ρ) ⊂ D met ρ := 1 −√ 1 − 2h
h α = 2α
1 +√
1 − 2h, dan heeft f een uniek nulpunt in
S(x0, r) ∩ D met r := 1 +√ 1 − 2h
n α
en convergeert het Newtonproces, startend in x0, naar dit nulpunt. (Geen bewijs).
Hoewel deze stelling op zich een mooi resultaat is, is de verficatie van (2.4a) in het algemeen ondoenlijk en kijkt men gewoon of het proces in de praktijk convergeert.
7.f Een aangepaste (gedempte) Newtonmethode
De Newtonmethode (7.16) convergeert lokaal, d.w.z. er is convergentie, als het startpunt voldoend dicht bij het gewenste nulpunt ligt. Starten we wat verder weg, dan kan de rij iteranden divergeren (of convergeren naar een andere limiet).
Voorbeeld: De funktie arctan(x) heeft precies een nulpunt, x = 0. De Newtoniteratie luidt:
xk+1 = xk− (1 + x2k) arctan (xk) . (7.18) Als voor het startpunt x0 geldt
arctan(x0) > 2 x0 1 + x20 ,
dan ligt de volgende iterand x1 verder van x = 0 dan dit startpunt, x1 = x0− (1 + x20) arctan (x0) < − x0.
Analoog geldt dan ∀k : |xk+1| > |xk|. Aangezien het proces (7.18) geen andere limiet kan hebben, moet de rij divergeren.
Een belangrijke reden, waarom (3.1) divergeert voor te grote startwaarden x0 is, dat de stap
−(1 + x2k) arctan(xk) weliswaar in de goede richting is, maar zijn doel voorbij schiet. Een kortere stap zou beter zijn: vervang (3.1) door
xk+1 = xk− λk(1 + x2k) arctan (xk) (7.19) en kies in iedere slag een dempingsfactor λk∈ [0, 1]. Als we λkgoed kiezen, zal de volgende iterand wel dichter bij het gezochte nulpunt komen. Als de iteranden eenmaal voldoend dichtbij zijn, wordt demping overbodig en kunnen we terugkeren naar het standaard (kwadratisch convergente) Newtonproces. Een gedempt n-dimensionaal Newtonproces zal dus luiden:
kies startwaarde x0; k := 0;
repeat yk:= f (xk);
Dk:= ∇f(xk); sk:= D−k1yk
kies dempingsparameter λk∈ (0, 1];
xk+1:= xk− λksk until convergentie
(7.20)
Aangezien we x willen vinden zo, dat f (x) = 0, lijkt het criterium “kf(xk)k voldoend klein” zeer geschikt, vooral als we de Euclidische norm nemen omdat
ϕ(xk) := kf(xk)k22 =
n
X
i=1
fi2 (xk) (7.21)
twee keer continu differentieerbaar is als f het is. Bovendien heeft deze functie de volgende prettige eigenschap ten opzichte van de Newton-zoekrichting sk in (7.20):
Stelling 7: Er is een λk> 0 zodat de funktie
g(t) := ϕ(xk− tsk) (7.22)
monotoon dalend is voor alle t ∈ [0, λk).
Bewijs: Het is voldoende te laten zien, dat g dalend is in t = 0, d.w.z. dat g′(0) < 0. We vinden g′(0) = d
dt ϕ (xk− tsk)|t=0 = d dt (
n
X
i=1
fi2(xk− tsk))| t = 0
= −
n
X
j=1
∂ϕ
∂xj
(xk) sk,j ( = −sTk ∇ϕ| x = xk)
= −
n
X
j=1 n
X
i=1
2 fi(xk) ∂fi
∂xj
(xk) sk,j
= −2 ykT Dk sk = − 2kykk2= −2g(0)
(7.23)
Het bepalen van een geschikte waarde van λk heet “lijnminimalisatie” (line search). In het eendimensionale voorbeeld zien we dat het vinden van een minimaliserende λk equivalent is met het oorspronkelijke nulpuntsprobleem. Het heeft dus geen zin om een exacte lijnminimalisatie te doen. Een eenvoudige strategie voor de keuze van λk wordt gegeven door de volgende algoritme:
λk := 1; r := ϕ(xk);
whileϕ(xk− λksk) ≥ r do λk:= 12λk end (7.24)
We beginnen dus bij de waarde λ = 1 (echte Newtonstap) en we halveren λ (geen afrondfout) totdat we in een punt komen met een “residu”, dat kleiner is dan ϕ(xk). Het is duidelijk, dat het zinvol is om met λ = 1 te beginnen, omdat dit op den duur (voor grote k) de optimale stapgrootte is.
Als we echter nog ver van ons doel verwijderd zijn, kan dit betekenen dat we in iedere slag van het gedempte Newtonproces (7.20) vele malen moeten halveren en ϕ(x − λs) berekenen bij iedere halvering. Dit werk kunnen we verminderen door in de k-de stap het zoeken naar een geschikte dempingsparameter niet aan te vangen met λ = 1, maar met tweemaal de demping van de vorige Newtonstap, tenzij dit groter dan 1 is:
if k = 0 then λ0 := 1 elsif λk−1 ≤ 1
2 then λk:= 2λk−1 else λk:= 1; end r := ϕ(xk);
while ϕ(xk− λksk) ≥ r do λk:= 12λk; end
(7.25)
Hierbij baseren we ons op de heuristiek dat “de dempingsparameter in de k-de slag wel niet veel van die van de k − 1-ste slag zal verschillen”, terwijl anderzijds de waarde van de dempingsparameter terug kan keren tot de optimale λ = 1 als we dicht genoeg bij het nulpunt zijn.
We kunnen in (7.24) en (7.25) beter gebruik maken van de informatie uit (7.23) over de afname van ϕ in de zoekrichting door niet slechts te eisen, dat ϕ(xk+1) < ϕ(xk) = g(0), maar te eisen dat een fractie α van de haalbare afname ook gerealiseerd wordt, dus
ϕ(xk+1) < g(0) + αλkg′(0) = g(0)(1 − 2αλk) met 0 < α < 1 . (7.26) In de praktijk blijkt α = 0.1 een goede keuze.
Als we g(1) = ϕ(xk − sk) berekend hebben maar de Newtonstap verworpen hebben omdat g(1) ≥ g(0)(1 − 2α), kunnen we echter een beter gebruik maken van de beschikbare gegevens. We kennen g(0), g′(0) en g(1) en het ligt dan voor de hand om de functie g te benaderen met de parabool
p(t) = g(0) − 2tg(0) + t2(g(1) + g(0)) , ¡waarbij p′(0) = g′(0) = −2g(0)¢, (7.27) welke in de genoemde punten met g overeenstemt; p is een kwadratisch “model” voor het verloop van ϕ langs de lijn xk− tsk. Deze parabool p(t) heeft altijd een minimum omdat de co¨effici¨ent van t2 positief is. Uit p′(t) = −2g(0) + 2t(g(0) + g(1)) = 0 volgt, dat dit minimum wordt aangenomen in het punt
t0:= g(0)
g(0) + g(1) ≤ 1 2(1 − α).
Als deze parabool een goed model is van g, zal de gevonden waarde van t0 een goede schatting leveren van het minimum van g en dus een betere waarde voor de dempingsparameter geven dan 12. Als we inderdaad g(t0) < g(0)(1 − 2αt0) vinden, hebben we een geschikte demping gevonden en kunnen we overgaan naar de volgende Newtonstap. Anders herhalen we deze procedure door naar het polynoom p1 te kijken, dat in t0 met g overeenstemt,
p1(t) := g(0) − 2tg(0) +t2
t20 (g(t1) + (2t0− 1) g(0)) . (7.28) Omdat g(t0) ≥ g(0)(1 − 2αt0), is de co¨effici¨ent van t2 positief en heeft p1 een minimum in
t1 := t20g(0)
g(t0) + (2t0− 1)g(0) ≤ t0 2(1 − α).
Als g(t1) < g(0)(1 − 2αt1) accepteren we de gevonden demping en anders herhalen we de procedure opnieuw.
Voor de afleiding van deze algoritme hebben we aangenomen, dat p een goed model is voor g.
Als echter g(1) ≫ g(0)) (of g(t0) ≫ g(0)), kun je hieraan twijfelen. In dat geval is de gevonden waarde van t0 veel kleiner dan 1 (resp. t1 ≪ t0). Het verdient dus aanbeveling om een drempel te leggen onder de te accepteren waarde van t0 (resp. t1), b.v. t0 > 0.1 (resp. t1 > 0.1 ∗ t0). Zo vinden we de algoritme:
tmin := max
µ g(0)
g(0) + g(1), 0.1
¶
; while g(tmin) ≥ g(0)(1 − 2αtmin) do
t1 := t2min g(0)
g(tmin) + (2tmin− 1) g(0) ; tmin:= max(t1, 0.1 ∗ tmin) end;
λk:= tmin
(7.29)
In deze algoritme wordt de staplengte steeds sterker verkleind dan bij de halveringsstrategie (7.24).
Omdat hierbij rekening gehouden wordt met de reeds gevonden waarden, mogen we verwachten, dat een goede demping zo sneller gevonden wordt.
7.g De methode van de steilste helling (steepest descent)
We beschouwen nu het probleem, een minimum te vinden van een functie ϕ(x) = ϕ(x1, . . . , xn) : IRn → IR, aangenomen dat zo’n minimum bestaat. De methode van Newton lost dit op door een nulpunt te bepalen van de gradient
∇ϕ(x) = µδϕ
δx1
, . . . , δϕ δxn
¶T
als funktie van x. Dit vereist kennis van de matrix van tweede afgeleiden van ϕ, de Hessiaan Hϕ := ∇∇Tϕ. Vaak is deze Hessiaan echter niet of slechts zeer moeilijk uit te rekenen, zodat we met een differentiebenadering moeten volstaan. We kunnen echter ook methoden zoeken, die geen gebruik maakt van de Hessiaan of benaderingen ervan.
We kunnen de funktie ϕ beschouwen als een “berglandschap” op een n-dimensionale aarde, waarin we moeten proberen het dal te vinden vanuit een zeker startpunt x0. We willen naar beneden lopen en kiezen hiervoor een richting s. Door de helling in deze richting te bepalen, kunnen we zien welke kant naar beneden gaat. We hebben het idee, dat we door steeds maar naar beneden te lopen uiteindelijk wel in het dal (het minimum) zullen aankomen. Dit is echter alleen waar als ons pad niet op den duur parallel gaat lopen met een horizontale richting, dus als de hoek met alle horizontale richtingen groter dan een zekere positieve waarde blijft. Het snelst (d.w.z. via de korste weg) zijn we beneden, als we steeds de richting van de steilste helling blijven volgen (met infinitesimaal kleine stapjes!). Deze steilste helling wordt gegeven door de gradient ∇ϕ van ϕ. De algoritme voor de gradient methode luidt dan:
kies startpunt x0; for k := 0, 1, 2, · · ·
sk := ∇ ϕ(xk); (bepaal de richting van de gradient)
bepaal λk zo dat ϕ(xk− λksk) ≤ ϕ(xk− tsk), ∀t , (lijnminimalisatie) xk+1:= xk− λksk
end
(7.30)
Hierbij dalen we dus steeds een eindje af in de richting van de steilste helling en wel zover, dat we het minimum in die richting bereiken. Dit proces blijkt te convergeren naar een stationair punt
van ϕ (punt met ∇ϕ = 0). Algemener zullen we bewijzen dat iedere afdaalrichting sk goed is, als de cosinus van de hoek tussen sk en de gradient ∇ϕ(xk) maar niet te klein wordt (d.w.z. als de hoek van de zoekrichting met een horizontale richting maar niet te klein wordt).
Stelling 8: Op een open gebied D ⊂ IRn is ϕ : D → IR een continu differentieerbare funktie, die voldoet aan de eis
k∇ϕ(x) − ∇ϕ(y)k2 ≤ γ kx − yk2, ∀x, y ∈ D . (7.31) Het punt x0∈ D is het beginpunt van het iteratieve proces
kies startpunt x0; for k := 0, 1, 2, · · ·
kies zoekrichting sk zo, dat voor zekere β > 0 geldt:
(a) sTk ∇ϕ(xk) ≥ β kskk2 k∇ϕ(xk)k2 ; kies stapgrootte λk zo, dat
(b) ϕ(xk− λksk) ≤ ϕ(xk− tsk), ∀t met xk− tsk ∈ D ; xk+1 := xk− λksk
end.
(7.32)
Als de verzameling K := {x ∈ D|ϕ (x) ≤ ϕ(x0)} een compact deel van D is, dan heeft de rij {xk} minstens ´e´en verdichtingspunt z en is ieder verdichtingspunt een stationair punt van ϕ.
Opmerkingen: Voorwaarde (7.32a) eist dat de cosinus van de hoek tussen de zoekrichting en de gradient niet te klein wordt; voorwaarde (7.32b) is de lijnminimalisatie die in iedere slag dient te worden uitgevoerd; eis (7.31) is niet nodig maar voor het bewijs wel gemakkelijk.
In (7.32b) wordt een exacte lijnminimalisatie gevraagd. In de praktijk is dit niet effici¨enten en doen we een “inexacte” lijnminimalisatie, waarbij we stoppen als een fractie van de haalbare vermindering van de objectfunctie bereikt is zoals in (7.26).
Bewijs: De rij {xk} is bevat in de compacte deelverzameling K en heeft dus minstens ´e´en verdicht- ingspunt z; we behoeven dus slechts te bewijzen, dat een verdichtingspunt een stationair punt is.
Bij definitie is de rij {ϕ(xk)} monotoon dalend; omdat ϕ continu is op een compacte verzameling K is de rij {ϕ(xk)} is naar beneden begrensd en dus convergent, d.w.z. limk→∞ ϕ (xk) = ϕ(z).
Stel nu dat z g´e´en stationair punt is, dan is er een α > 0 zodat k∇ϕ(z)k ≥ α > 0 en er is een bol S(z, r) met straal r rond z zodat
k∇ϕ(x)k ≥ 1
2α > 0 ∀x ∈ S(z, r) (7.33)
We kunnen nu een stukje Taylorontwikkeling bekijken van ϕ langs de zoekrichting xk−tsk. Hiervoor geldt wegens (7.31)
|ϕ(xk− tsk) − ϕ(xk) + tsTk ∇ϕ(xk)| ≤ 1 2 γ t2 als kskk2= 1, zodat
ϕ(xk− tsk) ≤ ϕ (xk) − tsTk ∇ ϕ(xk) +1 2 γ t2 Het rechterlid is minimaal als
t = tm:= 1
γ sTk ∇ϕ(xk) met minimum ϕ(xk) − 1
2γ (sTk ∇ϕ(xk))2.
In het volgende punt moet ϕ kleiner zijn dan dit minimum, dus ϕ(xk+1) ≤ ϕ(xk) − 1
2γ (sTk ∇ϕ(xk))2
≤ ϕ(xk) −β2
2γ k∇ ϕ(xk)k2 op grond van (7.32a)
≤ ϕ(xk) −α2β2
8γ op grond van (7.33).
Aangezien er oneindig veel punten van de rij {xk} binnen S(z, r) liggen, volgt hieruit dat de rij {ϕ(xk)} niet kan convergeren, hetgeen strijdig is.
Gevolg. Het gedempte Newtonproces (7.20) is convergent als λk wordt gekozen door lijnminima- lisatie; ga na dat de cosinus van de hoek tussen de zoekrichting en de gradient niet nul kan worden.
Opmerking: Bij ongunstige keuze van het startpunt kan het aantal iteraties bij een gradientmethode sterk oplopen, zoals in fig. 16 geschetst is. Getekend zijn de hoogtelijnen (niveaulijnen) behorende bij de iteranden xn, xn+1, xn+2· · · en de zoekrichtingen voor minimalisatie van de funktie f(x, y) :=
x2+ 9y2:
-5 -4 -3 -2 -1 0 1 2 3 4 5
0 2 4 6 8 10
gradient algorithme voor x^2 + 9 y^2 ; startpunt x=9, y=1
Figure 16: Voorbeeld van de iteratiestappen van de gradi¨entalgoritme toegepast op de funktie f (x, y) := x2+ 9y2 en startpunt (x0, y0) := (9, 1). Vanuit het startpunt wordt gezocht naar het minimum van f gezocht in de richting van de gradi¨ent, een lijn met helling 45o. In het minimum doen we een volgende stap in een richting loodrecht op de vorige. Zo zullen er dus zeer veel kleine stapjes nodig zijn om bij het minimum te komen.
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.
