Situation, in der beide Methoden extrem lang-sam konvergieren:
6.1.5 Newton-Raphson
Wenn die 1.Ableitung von f(x) analytisch zur Verf¨ugung steht, kann man statt lokaler Sekantenapproximationen auch die Tan-gente verwenden, um eine neue Approximation an die Nullstelle zu konstruieren. Da man f¨ur die Sekante zwei Punkte braucht, aber f¨ur die Tangente nur einen (und die Steigung), kommt man dann auch ganz ohne Intervallschachtelungen aus:
1. Rate einen Startpunkt xi, i= 1 2. Berechne f(xi) und f0(xi)
3. Approximiere die Funktion lokal durch die Tangente, d.h. die n¨achste Approximation f¨ur die Nullstelle ist die Nullstelle der Geraden durch den Punkt (xi, f(xi)) mit der Steigung f0(xi):
f(x) ≈ f(xi) +f0(xi)(x−xi) = 0! ⇒ xi+1 = xi − f(xi)
f0(xi) (126)
4. Wenn die Schrittweite abs(xi+1 −xi) oder der Funktionswert f(xi+1) kleiner ist als eine Toleranzgrenze, stop.
5. Sonst i ←i+ 1 6. Gehe zu (2).
Beachte: pro Iterationsschritt brauchen wir hier zus¨atzlich zur Funk-tionsberechnung auch einmal eine Ableitungsberechnung. Wenn diese Ableitung numerisch approximiert werden muß (finite Differenzen), ist es daher besser, regula falsi oder das Sekantenverfahren zu verwenden!
Vorteile:
• kein Anfangsintervall n¨otig, nur ein einzelner Startwert
• quadratische Konvergenz: i+1 = −2i2ff000(x)(x) Nachteile:
• die Tangentenextrapolation kann auch spektakul¨ar schiefgehen;
• keine M¨oglichkeit, die Nullstellenapproximationen in einem eingrenzenden Intervall zu halten.
Fehler im Newton-Verfahren: Oszillation:
Fehler im Newton-Verfahren: Divergenz:
bessere Strategie: mache Bisektion (oder Sekante/regula falsi), wenn die Newton-Schritte nicht schnell genug konvergieren oder aus dem vorgegebenen Intervall (bracket) hinausf¨uhren; sonst akzeptiere die Newton-Schritte.
6.2 Extremwertsuche = Optimierung, in 1D
Maximierung von f(x) ist Minimierung von −f(x), daher ist im folgenden immer nur von Minimierung die Rede.
Zur Eingrenzung von Minima braucht man nicht zwei, sondern drei Punkte:
Die stetige Funktion f(x) hat im Intervall [a, b] mit a < c < b ein Minimum, wenn f(c) < f(a) und f(c) < f(b).
Die drei Punkte a, b, c sind jedoch im Prinzip leichter zu finden als bei der Nullstellensuche: Die Vorschrift
”gehe runter (ggf.
in immer gr¨oßeren Schritten), bis die Funktion wieder ansteigt“, ist sehr oft erfolgreich, aber auch nicht immer: 2 Beispiele f¨ur Fehlersituationen:
1 2
3
4
1
2
3
4
6.2.1 Intervallschachtelung nach dem goldenen Schnitt
Minimierung in Analogie zur Bisektion bei der Nullstellensuche (ebenfalls mit einer Funktionsauswertung pro Iteration):
1. Gegeben sei ein Intervall [a, b] und ein weiterer Punkt c, die ein Minimum eingrenzen (s.o.).
2. Konstruiere einen neuen Punkt d, sodaß
• d im gr¨oßeren der beiden Intervalle [a, c] bzw. [c, b] liegt (hier: [c, b]),
• und zwar von c aus gesehen um die Strecke w×(b−c) in dieses Intervall hinein, mit w = (3−√ 5)/2
(dadurch ¨uber die sukzessiven Intervallschachtelungen skaleninvariant optimale Unterteilung, siehe Numerical Recipes) 3. Abbruch, wenn c, d ausreichend dicht beieinander; beste Approximation ans Minimum: der kleinere der beiden Werte
f(c), f(d). Sonst:
4. W¨ahle aus den jetzt vier Punkten a, c, d, b ein neues Triplett von Punkten, sodaß das neue Triplett das mit den niedrigst m¨oglichen Funktionswerten ist und gleichzeitig noch das Minimum eingrenzt.
5. Gehe zu (2).
Das jeweils n¨achste Intervall ist hier 0.61803 mal so breit wie das vorige; nicht ganz so gut wie das Verh¨altnis 0.5 bei Bisektion (in beiden F¨allen ist die Konvergenz linear).
6.2.2 Inverse parabolische Interpolation
Ahnlich wie eine Funktion in der N¨¨ ahe einer Nullstelle n¨aherungsweise linear ist, ist eine Funktion in der N¨ahe eines Minimums n¨aherungsweise parabolisch. Das Pendant zu den Sekantenmethoden ist daher die inverse parabolische Interpolation:
1. Gegeben sei ein Intervall [a, b] und ein weiterer Punkt c, die ein Minimum eingrenzen (s.o.).
2. Konstruiere eine Parabel durch a, c, b; der neue Punkt d ergibt das Minimum dieser Parabel und ist explizit gegeben durch:
d = b− 1 2
(b−a)2[f(b)−f(c)]−(b−c)2[f(b)−f(a)]
(b−a)[f(b)−f(c)]−(b−c)[f(b)−f(a)] (127) 3. Abbruch, wenn die Punkte c, d ausreichend dicht beieinander liegen; beste Approximation ans Minimum: der kleinere der
beiden Werte f(c), f(d). Sonst:
4. W¨ahle aus den jetzt vier Punkten a, c, d, b ein neues Triplett von Punkten, sodaß das neue Triplett das mit den niedrigst m¨oglichen Funktionswerten ist und gleichzeitig noch das Minimum eingrenzt.
5. Gehe zu (2).
Problematisch dabei ist, daß Gl.127 genausogut zu einem Maxi-mum wie zu einem MiniMaxi-mum f¨uhren kann und nat¨urlich versagt, wenn die Funktionsform stark von einer Parabel abweicht.
Ein robusteres Praxisverfahren ist daher die Minimierung nach Brent, bei der je nach Erfolg Schritte gem¨aß inverser paraboli-scher Interpolation bzw. gem¨aß goldenem Schnitt eingesetzt wer-den.
6.2.3 1D-Minimierung mit Gradienteninformation
Alle obigen 1D-Minimierungsverfahren ben¨otigen nur Funktionswerte, keine ersten (oder h¨oheren) Ableitungen. Wie bei der Nullstellensuche (Newton-Raphson), kann Ableitungsinformation n¨utzlich sein, hat aber auch Nachteile:
• ein einziger Ableitungswert entspricht im Informationsgehalt mehreren Funktionswerten,
• aber nur in einer hinreichend kleinen x-Umgebung;
• je nach Anwendungsfall kann das Berechnen von Ableitungen (im Vergleich zum Berechnen von Funktionswerten) einfach und billig oder aufwendig und teuer oder manchmal sogar unm¨oglich sein.
Daten aus ersten Ableitungen k¨onnen in unterschiedlicher Weise eingesetzt werden, zwischen folgenden Extremen:
• konservativ: nur zur Entscheidung, ob ein neuer Punkt d in [a, c] oder in [c, b] liegen sollte;
• aggressiv: Konstruktion von interpolierenden Polynomen h¨oherer Ordnung, unter Verwendung von Funktions- und Ablei-tungsinformation einiger alter Punkte.
Welche Strategierichtung erfolgreicher ist, h¨angt davon ab,
”wie pathologisch“ die untersuchte Funktion ist.
6.3 Extremwertsuche = Optimierung, in nD Anwendungsbeispiele:
• ”Geometrieoptimierung“ molekularer Systeme: optimiere potentielle Energie V durch Variation der Atompositionen 3 – Minima = (meta)stabile Molek¨ulkonfigurationen;
– Sattelpunkte 1. Ordnung = (m¨ogliche) ¨Ubergangszust¨ande chemischer Reaktionen.
• Hartree-Fock, DFT, CAS-SCF,. . . optimieren die elektronische Energie durch Variation der Wellenfunktion (der Orbitale) Gemeinsamkeiten von 1D- und nD-Minimierung:
• es gibt i.d.R. mehrere Minima; welches man findet, h¨angt v.a. vom Startpunkt ab,
• dies ist meist nur ein lokales Minimum, nicht das globale.
Unterschiede von nD-Minimierung gegen¨uber 1D-Minimierung:
• Anzahl Minima steigt i.d.R. exponentiell mit n.
• 1D: 2 Schrittrichtungen; bereits in 2D: ∞ viele Richtungen! Dadurch:
– keine strikte Eingrenzung m¨oglich.
– Gradienteninformation f¨ur Schrittrichtungswahl fast unverzichtbar
6.3.1 nD-Minimierung ohne Gradienten
• Neue Methoden: BOBYQA,4NEWUOA,5welche geschickt mehrdimensionale, quadratische Interpolation durchf¨uhren.
• Altere Algorithmen (Simplex, Powell, Brent): Siehe¨ Numerical Recipes.
3H. B. Schlegel,WIREs Comput Mol Sci, 1, 790-809 (2011).
4M. J. D. Powell, The BOBYQA algorithm for bound constrained optimization without derivatives, Cambridge NA Report NA200906, 2009
5M. J. D. Powell, Developments of NEWUOA for minimization without derivatives, IMA J. Numer. Anal., 28, 649–664 (2008).
6.3.2 Methode des steilsten Abstiegs (steepest descent)
1. Gegeben sei ein Startpunkt ~xi, i= 1, im Suchraum.
2. Berechne den Gradienten ~gi = ∇f(~xi); die Richtung des steilsten Abstiegs ist dann gegeben durch −~gi.
3. Verwende einen beliebigen 1D-Minimierungsalgorithmus, um in dieser Richtung das Minimum zu finden: minαf(~xi − α~gi) (line search, muß nicht sehr genau sein!).
4. Wenn |(f(~xi)−f(~xi−1))/f(~xi)| < Toleranz ϑ oder ||~gi|| < ϑ, stop.
5. Sonst: i ←i+ 1; gehe zu (2)
Vorteil: garantierter Fortschritt in jedem Iterationsschritt.
Nachteil: Sukzessive Schrittrichtungen stehen per Konstruktion aufeinander senkrecht. Das ist in sehr speziellen Situationen (sph¨arisch symmetrische Minimum-Umgebung, s.u.) gut, aber im Allgemeinen nicht:
• Schrittrichtungen zeigen nicht zum Minimum ⇒“Zickzack”-Wege;
• sp¨atere Schritte zerst¨oren teilweise den Optimierungsfortschritt in die fr¨uheren Richtungen.
6.3.3 nD-Minimierung via lokale quadratische Approximation
Approximation der Funktion f am Punkt ~x0 durch Taylorentwicklung bis zur zweiten Ordnung f(~x) ≈ f(~x0) +~g0T(~x−~x0) + 1
2(~x−~x0)TH0(~x−~x0). (128) (~g0: Gradient df /d~x am Punkt ~x0; H0: Hesse-Matrix d2f /d~x2 am Punkt ~x0)
Erste Ableitung von f:
∇f(~x) =~g0 +H0(~x−~x0). (129)
Nullsetzen der ersten Ableitung liefert:
~
x1 = ~x0 −H−10 ~g0 (130)
Dieser Newton-Schritt f¨uhrt von einem beliebigen Startpunkt~x0 ohne(!) line-search direkt ins Minimum der quadratischen Form.
Aber:
• H ist ggf. nicht (oder nicht exakt) verf¨ugbar, oder zu groß;
• H muß invertiert werden;
• die quadratische Form ist nur in einer kleinen Umgebung von ~x0 eine gute N¨aherung.
Zwei wesentliche Klassen von Verfahren: conjugate gradient, quasi-Newton.
Vorsicht in der Praxis:
• Newton-Schritte k¨onnen auch zu Sattelpunkten erster oder h¨oherer Ordnung f¨uhren (wg. numerischen Instabilit¨aten aber nur, wenn der Startpunkt nah genug am Sattelpunkt ist). Unterscheidung: Vorzeichen der Eigenwerte von H.
• Alle Methoden (auch steepest-descent) finden nur das Minimum, in dessen Einzugsbereich der Startpunkt liegt. Um andere Minima zu finden, muß man an anderen Stellen starten oder aus gefundenen Minima wieder entkommen.
conjugate gradient:
Steepest descent mit aufeinander senkrechten Schritten ~sTi+1~si = 0 w¨are gut (Schritte f¨uhren direkt zum Minimum und zerst¨oren Optimierungserfolg vorheriger Richtungen nicht) in einer quadratischen Form mit exakt gleichen Halbachsenl¨angen:
Im Allgemeinen hat eine quadratische N¨aherung jedoch eine ellipsoide Form und vormals zueinander senkrechte Vektoren sind dann zueinander konjugiert: ~sTi+1H~si = 0
Im conjugate-gradient-Verfahren gelingt es,
• Schritte in zueinander konjugierte Richtungen zu machen und dabei H indirekt zu ber¨ucksichtigen,
• ohne H jemals explizit berechnen oder speichern zu m¨ussen.
(siehe Numerical Recipes, Henrik Larssons Numerikskript WS15/16 sowie die sehr ausf¨uhrliche Erl¨auterung in: J. R. Shewchuk, An Introduction to the Conjugate Gradient Method without the Agonizing Pain, Edition 114, 1994,
https://www.cs.cmu.edu/∼quake-papers/painless-conjugate-gradient.pdf)
Daher ist conjugate-gradient z.B. gut f¨ur Strukturoptimierung sehr großer Molek¨ule, z.B. Proteine, bei denen die Berechnung der Hessematrix auf Kraftfeldniveau nicht extrem aufwendig ist, aber die Matrix schlicht zu groß zur Speicherung oder gar Inversion ist. Die Konvergenz zum Minimum erfolgt jedoch sehr viel schneller als bei steepest descent.
Ablaufschema von conjugate-gradient:
1. Der erste Schritt ist ein steepest-descent-Schritt: Berechne Gradientenvektor~g1 = ∇f~ (~x1) am Startpunkt~x1, dann line-search entlang ~s1 = −~g1, liefert neuen Punkt ~x2.
2. Berechne den aktuellen Gradientenvektor ~gn am aktuellen Punkt ~xn. 3. Berechne βn nach der Fletcher-Reeves-Formel:
βn = ~gnT~gn
~
gTn−1~gn−1 (131)
(bei steepest-descent sind alle βn = 0)
4. Erzeuge eine neue konjugierte Richtung ~sn = −(~gn +βn~sn−1) 5. line-search entlang ~sn, liefert neuen Punkt xn+1
6. gehe zu (2)
Die obige Formel f¨ur β ist die urspr¨unliche; andere, modifizierte Formeln sind f¨ur eine tats¨achliche quadratische Funktion ¨ aqui-valent, liefern aber bei einer Verwendung der quadratischen Form als Approximation i.d.R. bessere Resultate, z.B. die Polak-Ribi`ere-Formel:
βn = ~gTn(~gn −~gn−1)
~gTn−1~gn−1 (132)
quasi-Newton:
Schritte in konjugierte Richtungen, mit expliziter, sukzessiver Verbesserung (update) einergen¨aherten Hesseschen Matrix (z.B. am Anfang eine Einheitsmatrix).
• verschiedene update-Formeln in Verwendung, meistens BFGS (Broyden-Fletcher-Goldfarb-Shanno).
• ab und zu zu große Schritte ⇒ Beschr¨ankung auf vermuteten G¨ultigkeitsbereich der quadratischen N¨aherung (trust radius)
• teure line-searches k¨onnen durch viele Varianten mit intelligenter Schrittweitensteuerung ersetzt werden (augmented Hessian, rational function optimization, eigenvector following, geometry DIIS).
• ben¨otigt Speicher von der Ordnung O(Dim.2). Daher nicht verwendbar f¨ur hochdimensionale Probleme (Dim. > 1000)
• Am h¨aufigsten verwendet: low-storage BFGS (L-BFGS), welches die Speicherung von H durch Speicherung von vorherge-henden Vektoren der letzten ∼10 Schritte vermeidet.
6.4 Mehrdimensionale nichtlineare Nullstellensuche
Eine M¨oglichkeit, das Problem f~(~x) = 0 zu l¨osen, ist eine einfache Erweiterung von 1D-Newton-Raphson: Wir entwickeln f~ in eine Taylorreihe:
f~(~x+δ~x) =f~(~x) +Jδ~~x+O(δ~x2) (133) mit der Jacobi-Matrix J:
Jij = ∂fi
∂xj (134)
Wenn wir quadratische und h¨ohere Terme vernachl¨assigen (! lineare Approximation, i.Ggs. zur quadratischen bei der multidi-mensionalen Minimierung) und f~(~x+δ~x) =~0 setzen, finden wir als Korrekturschritt δ~x zur Nullstelle:
~
xi+1 = ~xi +δ~x = ~xi−J−1f~ (135)
Dummerweise ist die direkte Verwendung dieser Vorschrift in nD erheblich kritischer als die Verwendung von Newton-Raphson in 1D.
Die i.A. zu gef¨ahrlichen Newton-Raphson-Schritte lassen sich jedoch folgendermaßen entsch¨arfen: F¨ur die Hilfsfunktion F = 12f~·f~ ist der Schritt nach Gl.135 ein Schritt steilsten Abstiegs:
∇F ·δ~x = (f~J)·(−J−1f~) =−f~·f <~ 0 (136) Also kann man zun¨achst immer den Schritt nach Gl.135 versuchen. Wenn sich dadurch jedoch F nicht verkleinert, reicht es aus, in derselben Richtung weniger weit zu gehen (backtracking), bis sich eine akzeptable Schrittweite ergibt.
Warum ist nD-Nullstellensuche schwieriger als nD-Minimierung?
Wir k¨onnen hoffen,nUnbekannte mitnGleichungen bestimmen zu k¨onnen, aber es gibt keine Garantie, daß eine L¨osung existiert.
Bei der oben vorgestellten nD-Minimierung sind die Gradientenkomponenten voneinander abh¨angig. Das ganze Problem ist vorstellbar als die Suche nach einem tiefsten Punkt in einer h¨ugeligen Landschaft. Es gibt daher lokale Information (Gradient → steilster Abstieg), die zwangsl¨aufig zum Erfolg f¨uhren muß.
Im Gegensatz dazu sind n nichtlineare Gleichungen i.A. unabh¨angig voneinander. F¨ur nur zwei Funktionen f(x, y) = 0 und g(x, y) = 0 ist das wie unten abgebildet vorstellbar. L¨osungen sind an den Schnittpunkten der Null-Konturen f = 0 und g = 0, die v¨ollig unabh¨angig voneinander sind und jeweils auch aus mehreren unzusammenh¨angenden Teilen bestehen k¨onnen. Um festzustellen, ob es L¨osungen gibt und wenn ja, wie viele, m¨ußte man alle diese Teil-Konturen komplett nach Schnittpunkten absuchen. Die Situation wird dadurch noch weiter verkompliziert, daß die oben eingef¨uhrte Hilfsfunktion F multiple Minima hat, z.B. an allen Stellen, an denen sich f = 0 und g = 0 nahe kommen, wie etwa im Punkt
”M“ – dort befindet sich jedoch keine L¨osung des Gleichungssystems.