„Downhill-Simplex-Verfahren“ – Versionsunterschied
[gesichtete Version] | [gesichtete Version] |
TeesJ (Diskussion | Beiträge) typo, Lit +1 |
→Beispielimplementierung: Korrekturen, siehe [Diskussion Benutzer:PeterFrankfurt#Nelder-Mead] |
||
Zeile 137: | Zeile 137: | ||
innerhalb eines Parametersatzes: |
innerhalb eines Parametersatzes: |
||
Par0[NP+1][i]=(1 + NA) * Par0[NP+3][i] - NA * Par0[0][i]; |
Par0[NP+1][i]=(1 + NA) * Par0[NP+3][i] - NA * Par0[0][i]; |
||
fs= |
fs=FQS(NP+1); // und Fehler berechnen |
||
// Fallunterscheidung, ob Verbesserung erreicht |
// Fallunterscheidung, ob Verbesserung erreicht |
||
Zeile 146: | Zeile 146: | ||
| innerhalb eines Parametersatzes: |
| innerhalb eines Parametersatzes: |
||
| Par0[NP+2][i]=(1 + NC) * Par0[NP+1][i] - NC * Par0[NP+3][i]; |
| Par0[NP+2][i]=(1 + NC) * Par0[NP+1][i] - NC * Par0[NP+3][i]; |
||
| fs= |
| fs=FQS(NP+2); // und Fehler berechnen |
||
| |
| |
||
| wenn fs>=fmin // wieder Vergleich mit bisher bestem Punkt |
| wenn fs>=fmin // wieder Vergleich mit bisher bestem Punkt |
||
Zeile 179: | Zeile 179: | ||
| innerhalb eines Parametersatzes: |
| innerhalb eines Parametersatzes: |
||
| Par0[NP+2][i]= NB * Par0[0][i] + (1 - NB) * Par0[NP+3][i]; |
| Par0[NP+2][i]= NB * Par0[0][i] + (1 - NB) * Par0[NP+3][i]; |
||
| fs= |
| fs=FQS(NP+2); // und Fehler berechnen |
||
| // P** besser als bisher schlechtester Punkt? |
| // P** besser als bisher schlechtester Punkt? |
||
| wenn fs<fmax // besser als bisher schlechtester? |
| wenn fs<fmax // besser als bisher schlechtester? |
||
Zeile 190: | Zeile 190: | ||
| | innerhalb eines Parametersatzes j: |
| | innerhalb eines Parametersatzes j: |
||
| | Par0[j][i]=(Par0[j][i] + Par0[1][i]) / 2; |
| | Par0[j][i]=(Par0[j][i] + Par0[1][i]) / 2; |
||
\ \ fs= |
\ \ fs=FQS(j); // und Fehler berechnen |
||
// Schleifenende |
// Schleifenende |
||
Version vom 22. April 2010, 01:26 Uhr
Der Simplex-Algorithmus nach Nelder und Mead (Comp. J., vol. 7, 1965, p. 308) oder auch Downhill-Simplex-Verfahren oder manchmal auch einfach Simplex-Algorithmus ist im Unterschied zum Namensvetter für lineare Probleme (Simplex-Algorithmus) eine Methode zur Optimierung nichtlinearer Funktionen von mehreren Parametern. Er fällt in die Kategorie der Hillclimbing- oder Downhill-Suchverfahren. Angewendet werden kann er z. B. auch beim Kurvenfitten.
Grundlagen
![](https://upload.wikimedia.org/wikipedia/commons/thumb/0/09/Nelder_Mead1.gif/320px-Nelder_Mead1.gif)
![](https://upload.wikimedia.org/wikipedia/commons/thumb/9/96/Nelder_Mead2.gif/320px-Nelder_Mead2.gif)
Als Besonderheit kommt dieser Simplex-Algorithmus ohne Ableitungen der Funktion nach den Parametern aus, welche bei manchen Problemstellungen einfach nicht mit vertretbarem Aufwand berechenbar sind. Durch sinnvolle Vergleiche der Funktionswerte mehrerer Punkte im Parameterraum wird ähnlich wie bei einer Regula falsi mit Schrittweitensteuerung die Tendenz der Werte und der Gradient Richtung Optimum angenähert. Das Verfahren konvergiert nicht extrem schnell, ist aber dafür einfach und relativ robust.
Der Algorithmus basiert auf einem Simplex, dem einfachsten Volumen im N-dimensionalen Parameterraum, das aus N+1 Punkten aufgespannt wird. Im Eindimensionalen auf der Zahlengeraden ist das eine Strecke, im Zweidimensionalen ein Dreieck usw. Jeder Punkt entspricht dabei einem (Koordinaten-)Satz Parameter (bzw. Variablen), und zu jedem Punkt kann man einen Funktionswert berechnen, den man z. B. als Fehler- oder Kostenwert ansehen kann. Unter diesen N+1 Punkten ermittelt man dann den „schlechtesten“ und den „besten“. Von einem Iterationsschritt zum nächsten wird normalerweise der „schlechteste“ dieser Punkte durch einen neuen ersetzt, der hoffentlich „besser“ ist. Den „besten“ Punkt behält man als bisher beste Lösung immer bei.
Bei einer Visualisierung im Zweidimensionalen mit Höhenlinien für die zu optimierende Funktion sieht man (in den Bildern rechts) ein Dreieck, das sich zielstrebig dem Optimum nähert, indem es sich um eine seiner Seiten klappt, sich streckt oder in der Nähe des Optimums um dieses zusammenzieht. Siehe auch die Abbildung ganz unten.
Ohne Ableitungen entfallen diverse Fallen wie pathologisches Verhalten dieser Ableitungen bei Unstetigkeiten, so dass das Verfahren zwar im Allgemeinen nicht so schnell konvergiert (nämlich linear) wie Optimierungsverfahren mit Ableitungen, aber wesentlich robuster arbeitet. Die möglichen Fallen beim Simplex-Verfahren sind ähnlich wie bei den meisten Optimierungsverfahren unverhoffte lokale Nebenminima (bzw. -optima), auf die konvergiert werden kann, oder dass das Verfahren sich für einen oder mehrere Parameter zu früh auf einen konstanten Wert zusammenzieht.
Der Algorithmus selbst
Das Verfahren findet das Optimum einer Funktion mit N Variablen (Parametern) einer Funktion im N-dimensionalen Parameterraum. Ein Programm geht üblicherweise so vor, dass
(1) N+1 Anfangspunkte gewählt werden (Beispiel: man wählt einen Punkt als Startpunkt und ermittelt N weitere, indem man von diesem aus jeweils einen der Parameter um einen bestimmten Faktor − z. B. 1.1 − erhöht; alternativ: Man gibt für jeden Parameter zwei verschiedene Werte vor, berechnet den ersten Punkt mit allen Erstwerten und die anderen N Punkte auch mit den Erstwerten, nur dass jeweils einer der Parameter durch seinen Zweitwert ersetzt wird);
(2) die Funktionswerte an diesen Punkten berechnet werden (je nach Problemstellung ein Fehlerwert, Höhenwert, Kosten, Zeitaufwand usw.) und der beste und der schlechteste unter ihnen ermittelt wird;
(3) bei einem genügend guten Wert das Verfahren beendet wird, ansonsten (nur) den bisher schlechtesten der N+1 Punkte durch einen neuen ersetzen (wodurch ein neues Simplex-Volumen gebildet wird); und
(4) die Schleife mit Schritt (2) weitergeführt wird.
Das eigentliche Verfahren nach Nelder und Mead besteht in der Strategie, mit der in Schritt (3) der neue Simplex gebildet wird.
Nelder und Mead haben dazu folgende Regeln aufgestellt:
Immer den schlechtesten Punkt am (Mittelpunkt des) restlichen Simplex reflektieren („Reflexion“, unter Benutzung des Faktors alpha).
a) Wenn dieser Punkt besser ist als alle anderen, mache einen weiteren, noch größeren
Schritt (Faktor gamma) in dieselbe Richtung („Expansion“).
b) Wenn dieser Punkt immerhin eine Verbesserung bringt, ersetze den schlechtesten durch ihn und beginne den
Algorithmus wieder von vorn.
c) Wenn dieser Punkt eine weitere Verschlechterung bringt, berechne einen neuen Punkt, indem der schlechteste Punkt mit dem Faktor beta näher an den besten Punkt heran gerückt wird. Ist das erfolgreich, beginne den Algorithmus wieder von vorn („Kontraktion“).
d) Wenn dieser neue Punkt eine weitere Verschlechterung bringt, verwirf ihn und komprimiere
den ganzen Simplex um einen bestimmten Faktor beta um den bisher besten Punkt und
beginne den Algorithmus wieder von vorn („Komprimierung“).
Die Koordinatenveränderungen der Parameter während dieser Schritte werden mittels
der Nelder/Mead-Parameter alpha, beta und gamma berechnet, die
normalerweise als 1 („Reflexion“), 1/2 („Kontraktion“ und „Komprimierung“) und 2 („Expansion“) angesetzt werden.
Nach diesen Regeln wird die Iteration so lange weitergeführt, bis ein Konvergenzkriterium erfüllt ist. Der Simplex wird sich dabei in Richtung des Optimums verlagern und sich schließlich um dieses herum zusammenziehen.
Erweiterungen für Notfälle
Wie bei anderen Verfahren besteht immer die Gefahr, dass auf ein Nebenoptimum konvergiert wird. Diese Gefahr kann man durch mehrere Maßnahmen mindern:
- Nach einer bestimmten Anzahl von Schritten wird neu angesetzt. Dazu wird der bisher beste Punkt beibehalten und um ihn herum ein neuer Start-Simplex aufgebaut, indem bei jedem weiteren Punkt die Koordinate (nur) eines Parameters z. B. um 5 % variiert wird.
- Genau wie eben, nur werden die zusätzlichen Simplex-Punkte durch Zufallswerte noch mehr flexibilisiert.
- Statt solche Maßnahmen nach einer festen Anzahl von Schritten zu ergreifen, kann man den Zeitpunkt auch vom aktuellen Stand der Iteration abhängig machen. Dazu muss das Problem natürlich besser bekannt sein, dass man sich dessen auch sicher sein kann.
- Alternativ kann man sogar den Zeitpunkt der Neuansetzung per Zufallszahl variabel gestalten.
- Eine zusätzliche Sicherheitsabfrage könnte den Fall abfangen, dass sich nur einer der Parameter vorzeitig festgelaufen hat, dann könnte man ihm eine Spezialbehandlung zukommen lassen.
All diese Varianten erfordern eine intensive Beschäftigung mit dem jeweiligen konkreten Problem, es gibt leider keine allgemein gültigen Kochrezepte dafür.
Erweiterung zum Kurvenfitten
Will man nicht eine einzelne analytische Funktion von N Parametern (Koordinaten) optimieren, sondern eine Theoriefunktion an eine Messkurve anpassen, ist nur ein kleiner zusätzlicher Schritt notwendig: Als zu optimierende Funktion verwendet man die Fehlerquadratsumme, also die Summe aus den Quadraten der Differenzen aus den einzelnen Messwerten (Kurvenpunkten) und den Theoriefunktionswerten an denselben Stellen, wobei letztere einerseits von der x-Koordinate der Messkurve abhängen und zusätzlich von den N Parametern. Ob man aus der Fehlerquadratsumme noch die Wurzel zieht, ist für das Simplex-Verfahren irrelevant.
Beispielimplementierung
// Deklarationen bool fertig=false; // Schleifenabbruch // Benutzervorgaben, -eingaben int NP= ... ; // Anzahl der Parameter, die in FQS() eingehen, // auch "Dimension des Parameterraums" double ParP[NP]={ ... }; // primärer Parametersatz double ParS[NP]={ ... }; // sekundärer Parametersatz, gibt Schwankungsbreite vor double FQS(int pf) // Zu minimierende Funktion, hängt { // von den Parametern Par0[pf][...] ab. ... // Muss je nach Problemstellung erstellt werden. // So schreiben, dass sie für den Parametersatz Nr. pf } // berechnet und am Schluss in Par0[pf][0] (s. u.) // abgespeichert wird. // Ende Benutzervorgaben // Initialisierung // Simplex-Parameter double NA = 1, NB = .5, NC = 2; // Alpha, Beta, Gamma // Parametersätze, NP+1 Stück für den Simplex selbst, plus 3 weitere double Par0[NP+3][NP]; // darin 1. Index: 0..NP die (NP+1) Parametersätze des Simplex, // dabei auf 0 immer der schlechteste Punkt des Simplex, // auf 1 immer der beste Punkt des Simplex // NP+1 Punkt P*, s. u. // NP+2 Punkt P**, s. u. // NP+3 Mittelpunkt P-Quer des Simplex, s. u. // darin 2. Index: 0 Fehlerwert für diesen Satz // 1..NP einzelne Parameterwerte dieses Satzes // Initialisierung // Simplex aus Primär- und Sekundärwerten Parametersätze 0 bis NP mit Primärwerten belegen; In Parametersätzen 1 bis NP (i) jeweils Parameter Nr. i mit Sekundärwert belegen; // Berechnung der Anfangs-Fehlersummen für Parametersätze 0 bis NP: feh=FQS(i); // speichert Ergebnis in Par0[i][0], s. o. // Schleife für Iteration, Hauptschleife while (!fertig) NT++; // Simplex-Schritt Nr. NT // schlechtesten und besten Punkt ermitteln (höchster bzw. niedrigster Fehlerwert) Werte fmin und fmax auf Fehlerwert von Satz Nr. 0 vorbelegen; für alle Parametersätze von 1 bis NP: wenn Fehler für Satz Nr. i größer als fmax: noch schlechterer Punkt, fmax auf diesen Wert, Punkt mit Nr. 0 tauschen wenn Fehler für Satz Nr. i kleiner als fmin: noch besserer Punkt, fmin auf diesen Wert, Punkt mit Nr. 1 tauschen // Konvergenzabfrage, sind wir fertig? // Dies muss auch je nach Problem individuell erstellt werden, // hier als Beispiel die Bedingung, dass alle Parameter nur noch um 1 Promille // schwanken. fertig=true; für alle Parameter in allen Parametersätzen von 0 bis NP: wenn nur für einen Parameter die relative Abweichung des Wertes zwischen schlechtestem und bestem Punkt größer als 0.001, dann: fertig=false; // doch noch nicht fertig // Mittelpunkt P-QUER des Simplex nach Index NP+3 Mittelwerte über Sätze 1 bis NP (also ohne schlechtesten Punkt Nr. 0!) bilden und als Satz Nr. NP+3 speichern; // Punkt P* durch Reflexion des schlechtesten Punkts // am Mittelpunkt P-QUER (mit Nelder/Mead-Parameter alpha=NA) nach Index NP+1 innerhalb eines Parametersatzes: Par0[NP+1][i]=(1 + NA) * Par0[NP+3][i] - NA * Par0[0][i]; fs=FQS(NP+1); // und Fehler berechnen // Fallunterscheidung, ob Verbesserung erreicht wenn fs<fmin // neuer Punkt P* im Vergleich zu bisher bestem / // ja, besser! | // weiteren Schritt zu Punkt P** in dieselbe Richtung nach Index NP+2 | // (mit Nelder/Mead-Parameter gamma=NC) | innerhalb eines Parametersatzes: | Par0[NP+2][i]=(1 + NC) * Par0[NP+1][i] - NC * Par0[NP+3][i]; | fs=FQS(NP+2); // und Fehler berechnen | | wenn fs>=fmin // wieder Vergleich mit bisher bestem Punkt | / // keine weitere Verbesserung: P** verwerfen, P* nach Index 0, | | // also bisher schlechtesten Wert durch neuen (etwas besseren) ersetzen | \ Parametersatz Nr. NP+2 nach Nr. 0 kopieren; | else // von fs>=fmin | / // ja, auch eine Verbesserung! | | // auch besser als P* ? | | wenn fs>Par0[NP+1][0] | | / // nein, P* weiterverwenden und damit bisher schlechtesten auf Index 0 ersetzen | | \ Parametersatz Nr. NP+1 nach Nr. 0 kopieren; | | else // von fs>Par0[NP+1][0] | | / // ja, P** weiterverwenden und damit bisher schlechtesten auf Index 0 ersetzen \ \ \ Parametersatz Nr. NP+2 nach Nr. 0 kopieren; else // von fs<fmin / // nicht besser als bisher bester, nun prüfen, ob P* wenigstens überhaupt eine Verbesserung | wenn fs von P* für wenigstens einen der Punkte 0..NP kleiner ist, dann: | / // ja, besser als mindestens einer, mit P* den bisher schlechtesten Punkt auf Index 0 ersetzen | \ Parametersatz Nr. NP+1 nach Nr. 0 kopieren; | else // von fs von P* | / // nein, weiterhin schlechter als die anderen Punkte | | // Zusatzprüfung: sogar schlechter als bisher schlechtester Punkt fmax? | | wenn fs>fmax // schlechter als bisher schlechtester? | | // ja, erstmal nichts tun | | else // von fs>fmax | | / // nein, also wenigstens bisher schlechtesten Punkt damit ersetzen | \ \ Parametersatz Nr. NP+1 nach Nr. 0 kopieren; | // neuer Punkt P** durch Kontraktion zwischen bisher schlechtestem Punkt (Index 0) | // und P-QUER (Index NP+3) nach Index NP+2 (mit Nelder/Mead-Parameter beta=NB) | innerhalb eines Parametersatzes: | Par0[NP+2][i]= NB * Par0[0][i] + (1 - NB) * Par0[NP+3][i]; | fs=FQS(NP+2); // und Fehler berechnen | // P** besser als bisher schlechtester Punkt? | wenn fs<fmax // besser als bisher schlechtester? | / // ja, bisher schlechtesten durch P** ersetzen | \ Parametersatz Nr. NP+2 nach Nr. 0 kopieren; | else // von fs<fmax | / // nein, keinerlei Verbesserung, Notfall | | // folgt allgemeine Kompression um Faktor 2 um den bisher besten Punkt herum | | für alle Parametersätze außer Nr. 1 (bisher bester): // also j=0 und 2..NP | | innerhalb eines Parametersatzes j: | | Par0[j][i]=(Par0[j][i] + Par0[1][i]) / 2; \ \ fs=FQS(j); // und Fehler berechnen // Schleifenende
Eine direkt lauffähige Implementierung in QBasic ist unter [1] zu finden.
Grafische Darstellung eines Beispiellaufs
Das Bild zeigt das Ergebnis eines Programmlaufs. Die roten und blauen Linien stellen die Gefällelinien für Täler bzw. Berge dar. Es handelt sich um einen Kreisgraben in einer schrägen Ebene, es interessiert hier das lokale Minimum innerhalb des Kreisgrabens. Die Simplex-Iteration findet das Minimum an der Kreuzungsstelle der roten Linien. In x- und y-Richtung reicht die Grafik von 0 bis 4.
Literatur
- Nelder, J. A.; Mead, R. A: A Simplex Method for Function Minimization. Computer Journal, 1965, 7, 308-313
- Lagarias, J. C.; Reeds, J. A.; Wright, M. H.; Wright, P. E.: Properties of the Nelder-Mead Simplex Method in low Dimensions. SIAM Journal on Optimization, 1998, 9, 112-147 [1]
- Press, W. H.; Teukolsky, S. A.; Vetterling, W. T.; Flannery, B. P.: Numerical Recipes in C. Cambridge University Press 1997
- Larry Nazareth, Paul Tseng: Gilding the Lily: A variant of the Nelder-Mead algorithm based on Golden-Section Search. Computational Optimization and Applications, Bd. 22 (2002), Nr. 1, S. 133–144
- Saša Singer, Sanja Singer: Efficient implementation of the Nelder-Mead search algorithm. Applied Numerical Analysis and Computational Mathematics, Bd. 1 (2004), Nr. 2, S. 524–534 (Zusammenfassung)
- Marco A. Luersen, Rodolphe Le Riche: Globalized Nelder–Mead method for engineering optimization. Computers and Structures 82 (2004), S. 2251–2260 (pdf-Datei)
- Andrej Prša, Tomaž Zwitter: Introducing adapted Nelder & Mead's downhill simplex method to a fully automated analysis of eclipsing binaries. In: C. Turon, K. S. O'Flaherty, M. A. C. Perryman (Hrsg.): The Three-Dimensional Universe with Gaia (ESA SP-576). Proc. Gaia Sympos. 2004, ESA Publications Division, Nordwijk 2005, S. 611–614
- David Byatt, Ian Coope, Chris Price: 40 years of the Nelder-Mead algorithm (pdf-Datei)