Matroids Matheplanet Forum Index
Moderiert von matroid
Mathematik » Numerik & Optimierung » Optimierung in 2 Variablen - Numerische Probleme oder nicht?
Autor
Kein bestimmter Bereich J Optimierung in 2 Variablen - Numerische Probleme oder nicht?
Delastelle
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 17.11.2006
Mitteilungen: 2209
  Themenstart: 2022-10-24

Hallo Leute! Wer möchte kann mal folgende Funktion optimieren (minimieren): ff = max( max(x1^2+x2^4,(2-x1)^2+(2-x2)^2),2*exp(-x1+x2) ); // ein Maximum aus 3 Funktionen bei Verfahren kann der Startpunkt x1 = 1, x2 = -0.1 benutzt werden. (Ich hätte es auch als Knobelaufgabe stellen können...) Viele Grüße Ronald


   Profil
hyperG
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 03.02.2017
Mitteilungen: 1923
  Beitrag No.1, eingetragen 2022-10-24

Hi, wie man am 3D-Plot schön sehen kann, wird der größte Teil der Funktion durch x1^2 + x2^4 dominiert: https://matheplanet.com/matheplanet/nuke/html/uploads/b/47407_max_Plot.png Die exp-Funktion trägt erst mit x2 > 9 und sehr kleinem x1 zur Max-Bildung bei. Wenn man hier z.B. 2 Variablen mit über 1 Mio. Nachkommastellen genau berechnen will, bringt es zeitlichen Vorteil, wenn man zunächst mit "einfacher Genauigkeit" {double oder auch if ...} abfragt, in welchem Bereich man ist, statt 3 mal 3 Ergebnisse auf über 1 Mio. Stellen zu berechnen. Das (2 - a)^2 + (2 - b)^2 = 8+(a-4)*a + (b-4)*b bringt numerisch nicht viel. Aber ich vermute, dass Du nur symbolisch optimieren willst. Da habe ich nichts gefunden. Es gibt zwar eine Wandlung der Max-Funktion in die Vorzeichenfunktion {sgn() }, aber das wird zu komplex: 1/4 sgn(b^4 - (2 - a)^2 + a^2 - (2 - b)^2) b^4 + 1/4 sgn(b^4 - (2 - a)^2 + ... {über 1000 Zeichen!} Zusatz: Oder wolltest Du die Ähnlichkeit der letzten beiden Funktionen gegenüber einer Geraden für große Argumente optimieren? https://matheplanet.com/matheplanet/nuke/html/uploads/b/47407_VergleichzurGerade.png Grüße


   Profil
Delastelle
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 17.11.2006
Mitteilungen: 2209
  Beitrag No.2, vom Themenstarter, eingetragen 2022-10-24

Hallo hyperG! Es ist angedacht, die Funktion numerisch zu lösen. Das Subgradientenverfahren soll Probleme in der Nähe der Lösung haben. Leider habe ich den Code zum Subgradientenverfahren - was ich früher mal gerechnet habe - nicht mehr gefunden. Ein SQP Verfahren findet: \sourceon Octave X = 1.1391 0.8995 OBJ = 1.9522 INFO = 104 ITER = 14 NF = 96 LAMBDA = [](0x1) \sourceoff Info 104 bedeutet zuletzt zu kleine Schrittweite beim SQP verfahren. Das Beispiel war aus einem alten Skript u.a. mit konvexer Optimierung und nicht überall differenzierbarer Zielfunktion. Viele Grüße Ronald


   Profil
hyperG
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 03.02.2017
Mitteilungen: 1923
  Beitrag No.3, eingetragen 2022-10-24

Ich vermute einen Schreibfehler in der Funktion, denn ich sehe da kein lokales Maximum! https://matheplanet.com/matheplanet/nuke/html/uploads/b/47407_3D_ExpFalsch.png Vermutung: beim exp() scheint entweder Vorzeichen falsch exp(-x1-x2) oder fehlt Klammer exp(-(x1+x2)) damit es logisch wird? Oder andere Definition für max()?


   Profil
Delastelle
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 17.11.2006
Mitteilungen: 2209
  Beitrag No.4, vom Themenstarter, eingetragen 2022-10-24

Hallo hyperG! Wie ist bei Dir der Funktionswert für das beim SQP-Verfahren gefundene Minimum? f(1.1391,0.8995) = 1.9522 Ansonsten muss ich noch mal die Funktion mit meinen Aufzeichnungen vergleichen... Bei Octave/contourf erhalte ich: https://www.matheplanet.de/matheplanet/nuke/html/uploads/b/15578_Bundle_Contourf.png Viele Grüße Ronald Edit: \showon \sourceon Octave MAXITER = 100; x0 = [1 -0.1]; [X, OBJ, INFO, ITER, NF, LAMBDA] = sqp(x0,'bundle',[],[],[],[],MAXITER) *** function [ff] = bundle(x) x1 = x(1); x2 = x(2); ff = max( max(x1^2+x2^4,(2-x1)^2+(2-x2)^2),2*exp(-x1+x2) ); % x* = (1.1390, 0.8996) % f(x*) = 1.952 % Startpunkt 1,-0.1 *** % plot bundle merk = zeros(1000,3); nr = 0; for x = 0:0.1:2 for y = 0:0.1:2 nr = nr + 1; ff = bundle([x y]); merk(nr,1:3) = [x y ff]; end end nr hold off plot3(merk(1:nr,1),merk(1:nr,2),merk(1:nr,3),'o'); xx = 0:0.1:2; yy = 0:0.1:2; [X,Y] = meshgrid(xx); Z = max( max(X.^2+Y.^4,(2-X).^2+(2-Y).^2),2.*exp(-X+Y) ); figure contourf(X,Y,Z); \sourceoff \showoff


   Profil
zippy
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 24.10.2018
Mitteilungen: 4266
  Beitrag No.5, eingetragen 2022-10-24

Ihr redet aneinander vorbei: \quoteon(2022-10-24 20:42 - hyperG in Beitrag No. 3) Maximum \quoteoff \quoteon(2022-10-24 20:52 - Delastelle in Beitrag No. 4) Minimum \quoteoff --zippy


   Profil
hyperG
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 03.02.2017
Mitteilungen: 1923
  Beitrag No.6, eingetragen 2022-10-25

\quoteon(2022-10-24 21:50 - zippy in Beitrag No. 5) Ihr redet aneinander vorbei: \quoteon(2022-10-24 20:42 - hyperG in Beitrag No. 3) Maximum \quoteoff \quoteon(2022-10-24 20:52 - Delastelle in Beitrag No. 4) Minimum \quoteoff --zippy \quoteoff Danke Zippi bei der Hilfe zur Findung der gemeinsamen Sprache :-) Numerische Methoden unter Mathematica: "NelderMead" https://en.wikipedia.org/wiki/Nelder%E2%80%93Mead_method {f(x1,x2)=1.95222449387091,{x1->1.139037550752,x2->0.899560017603963}} MeshSearch = 1.9522244938710762 RandomSearch = 1.9522244985340693 die anderen 3 sind noch schlechter bei voreingestellter Genauigkeit & MaxIteration. (lassen sich also verbessern)


   Profil
Delastelle
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 17.11.2006
Mitteilungen: 2209
  Beitrag No.7, vom Themenstarter, eingetragen 2022-10-25

Hallo, ich weiß noch, dass ein Subradientenverfahren nach 700.000 Iterationen noch nicht am Ziel war. Das Verfahren konnte immer nur winzige Schritte in eine Abstiegsrichtung machen. Aber mein damaliges Programm habe ich wohl nicht aufgehoben. Vielleicht bekomme ich es wieder hin! Dazu kann man ja einen exakten Gradienten für die 3 Teilfunktionen angeben. Die aktuelle Funktion soll ein Beispiel für Probleme eines Abstiegsverfahrens (Subgradienten-Verfahren) sein. Nelder-Mead verwendet ja keine Ableitungen und kann bei "einfachen" Problemen durchaus schnell sein. Überrascht hatte mich, dass unter Octave das eingebaute SQP (Sequentiell Quadratic Programming) keine Probleme bei dem Startwert von (1,-0.1) hat. Da hatte ich mehr Schwierigkeiten erwartet. Eventuell kann ich ja mal ein einfaches Gradienten-Verfahren mit numerischen Gradienten und einfachen Schrittweiten anwenden. Mal sehen wie gut das geht. Oder ob das fast analog zum früheren Subgradientenverfahren wäre... Viele Grüße Ronald


   Profil
hyperG
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 03.02.2017
Mitteilungen: 1923
  Beitrag No.8, eingetragen 2022-10-25

\quoteon(2022-10-25 02:03 - Delastelle in Beitrag No. 7) Hallo, ich weiß noch, dass ein Subradientenverfahren nach 700.000 Iterationen noch nicht am Ziel war... \quoteoff Die Definition für "AM ZIEL" ist hier sehr wichtig! Genau genommen müsste man zunächst die Randbedingungen exakt vorgeben: -Suchbereich - Nachkommastellen-Genauigkeit beim Endergebnis (nicht bei den Zwischenergebnissen) - genaue Definition, ob lokales oder globales Min/Maximum - einige Algorithmen benötigen auch definierten Startwert Die 4 Nachkommastellen sind aber wirklich wenig!? Rechnet Octave nur mit single-Genauigkeit (float)? Da ist man ja mit primitivem "Probieren" oder "Random" schnell am Ziel... Oberhalb 16 NK brauchen einige Algorithmen ewig lange. Hier mal 2 Beispiele für gute & schlechte Parameter von Numerischer Min-Suche: \sourceon mathematica FindMinimum...Method -> "InteriorPoint", MaxIterations -> 44000, WorkingPrecision -> MachinePrecision (*, also etwa 15 Nachkommastellen braucht mit Startwert 1.14, 0.9 ()also schon relativ dicht am Ziel über 148 s und ist beim Endergebnis mit 1.9524... gerade mal 3 Nachkommastelle richtig -> sehr schlechter Algorithmus hingegen *) NMinimize ... Method -> "NelderMead" hat nach 0.03 s schon 13 richtige Nachkommastellen \sourceoff Übrigens fand ich auch interessant, dass man mit Hilfe von Mathematica aus a^2 + b^4 == (2 - a)^2 + (2 - b)^2 nach Umstellung b= eine Gleichung 4. Grades umgestellt nach Wurzelformel aus über 700 Zeichen machen konnte. So müsste auch die 1. MAX-Funktion eliminierbar sein... (bringt aber Rechentechnisch mehr Rundungsfehler ins System) Interessant: der 3. Befehl Minimize {also eigentlich Nicht Numerisch} \sourceon nameDerSprache 50 Stellen genaue Variablen: 1.952224493870660... mit Exp-Anteil nur 13 Stellen genaues Endergebnis 1.952224493870658... ohne Exp-Anteil -> scheint wichtig zu sein, oder? 100 Stellen genaue Zwischenergebnisse: 1.9522244938706589943... mit Exp-Anteil nur 15 Stellen genaues Endergebnis 1.9522244938706589939... 200 Stellen: WorkingPrecision -> 200 1.9522244938706590... wieder größer -> da stimmt was nicht!!! ah, doch: N[Minimize[\[Ellipsis]]] calls NMinimize for optimization problems that cannot be solved symbolically. \sourceoff


   Profil
hyperG
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 03.02.2017
Mitteilungen: 1923
  Beitrag No.9, eingetragen 2022-10-26

Da mich die große Ungenauigkeit bei Mathematicas numerischer Minimierung störte, habe ich mal mit Tricks symbolisch versucht und landete beim neuen Minimum: 1.952224493870658993966607731144527397669385249689703872043955198407665037355128618... was man in wenigen ms schnell berechnen kann! Werde ich mal genauer überprüfen...


   Profil
Delastelle
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 17.11.2006
Mitteilungen: 2209
  Beitrag No.10, vom Themenstarter, eingetragen 2022-10-26

Hallo, das Subgradientenverfahren was ich oben erwähnte, sieht so aus: https://www.matheplanet.de/matheplanet/nuke/html/uploads/b/15578_Subgradientenverfahren.jpg Eine größere Genauigkeit als ein paar Stellen strebe ich momentan nicht an. Bei vielen Problemen ist man zufrieden, wenn man ein paar Stellen berechnet hat. Viele Grüße Ronald Edit: ich habe die Funktion unter Fortran mit einem einfachen Gradientenverfahren mit numerischen Gradienten getestet. Dabei kam ich nach ca. 130.000 Iterationen auf eine relativ gute Lösung. Die Schrittweite zum Gradienten wurde sehr einfach ermittelt: aus (1, 0.5, 0.25, 0.125, ... ) // 40 Stück wird die Schrittweite genommen für die der Zielfunktionswert am niedrigsten ist. Iter = 130000 x1 = 1.18658178749 x2 = 0.861297350337 ff = 1.95829291280 Zeit = 0.358802 \showon \sourceon Fortran program grad c x* = (1.1390, 0.8996) c f(x*) = 1.952 implicit none double precision feld(1000000) integer lala,maxiter,i,og,ausgabeiter double precision x1,x2,x1best,x2best,x1aktuell,x2aktuell double precision x1test,x2test,kkkgenommen double precision x1plus,x2plus,ff,ff1,ff2,ffmin,h double precision grad99,grad1,grad2,maxstep,kkk real t0,t1 call cpu_time(t0) maxiter = 150000 ausgabeiter = 10000 x1 = 1.0d0 x2 = -0.1d0 og = 40 x1test = 1.1390d0 x2test = 0.8996d0 call bundle(x1test,x2test,ff) print *,'ff(opt) = ',ff do 100 lala = 1,maxiter call bundle(x1,x2,grad99); h = 0.0000001d0 x1plus = x1+h call bundle(x1plus,x2,grad1) ff1 = (grad1-grad99)/h c print *,ff1 x2plus = x2+h call bundle(x1,x2plus,grad2); ff2 = (grad2-grad99)/h c print *,ff2 maxstep = 2.0d0; do 200 i = 1,og kkk = maxstep/(2.0d0**i); c print *,'step ',kkk x1aktuell = x1 - kkk*ff1; x2aktuell = x2 - kkk*ff2; if (i == 1) then call bundle(x1aktuell,x2aktuell,ffmin); x1best = x1aktuell x2best = x2aktuell kkkgenommen = kkk else call bundle(x1aktuell,x2aktuell,ff); if (ff < ffmin) then ffmin = ff x1best = x1aktuell x2best = x2aktuell kkkgenommen = kkk endif endif 200 continue x1 = x1best x2 = x2best feld(lala) = kkkgenommen c print *,lala,kkkgenommen if (mod(lala,ausgabeiter) == 0) then call cpu_time(t1) print *,' ' print *,'Iter = ',lala print *,'x1 = ',x1 print *,'x2 = ',x2 print *,'ff = ',ffmin print *,'Zeit = ',t1 endif 100 continue c call ausgabe(feld,maxiter) call cpu_time(t1) print *,'Zeit ',t1 end c *********** c c *********** subroutine bundle(x1,x2,ff) implicit none double precision ff1,ff2,ff3,x1,x2,ff ff1 = x1**2.0d0+x2**4.0d0 ff2 = (2.0d0-x1)**2+(2.0d0-x2)**2.0d0 ff3 = 2.0d0*exp(-x1+x2) ff = max(ff1,ff2) ff = max(ff,ff3) end c *********** c c *********** subroutine ausgabe(feld,maxiter) implicit none double precision feld(1000000) integer i,maxiter c open(11,file='kurven.dat') do 100 i = 1,maxiter write(11,*) feld(i) 100 continue close(11) return end \sourceoff \showoff


   Profil
Delastelle
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 17.11.2006
Mitteilungen: 2209
  Beitrag No.11, vom Themenstarter, eingetragen 2022-10-26

Ich habe jetzt ein Subgradientenverfahren. Bei mir braucht es für das Finden einer Stelle nahe dem Optimum bis zu 1 Mio. Iterationen. Dazu wird der exakte Gradient verwendet. In Abhängigkeit welche der 3 Funktionen aktiv ist, wird Gradient 1 bis 3 verwendet. Es funktioniert einigermaßen für dieses Beispiel. Bei og = 15 (größere Schrittweite) auch "schon" nach ca. 30000 Iterationen. (Ausgabe bei 1 Mio. Iterationen) Iter = 1000000 x1 = 1.13951967768 x2 = 0.899182528422 ff = 1.95222549082 Schrittweite = 1.907348632813E-06 Zeit = 1.27921 \showon \sourceon Fortran program subgrad c x* = (1.1390, 0.8996) c f(x*) = 1.952 implicit none double precision feld(1000000) integer lala,maxiter,i,og,ausgabeiter,zustand,nr double precision x1,x2,x1best,x2best,x1aktuell,x2aktuell double precision x1test,x2test,kkkgenommen double precision ff,ff1,ff2,ffmin double precision maxstep,kkk,hurz real t0,t1 call cpu_time(t0) maxiter = 1000000 ausgabeiter = 50000 x1 = 1.0d0 x2 = -0.1d0 og = 20 c og = 25 // 50 Mio ff = 1.95 x1test = 1.1390d0 x2test = 0.8996d0 call bundle(x1test,x2test,ff,zustand) print *,'ff(opt) = ',ff do 100 lala = 1,maxiter call bundle(x1,x2,ff,zustand) c print *,zustand c h = 0.0000001d0 c x1plus = x1+h c call bundle(x1plus,x2,grad1) c call random_number(hurz) c zustand = floor(3*hurz)+1 c print *,zustand nr = 1 call subgrad99(x1,x2,ff1,zustand,nr) c ff1 = (grad1-grad99)/h c print *,ff1 c x2plus = x2+h c call bundle(x1,x2plus,grad2); nr = 2 call subgrad99(x1,x2,ff2,zustand,nr) c ff2 = (grad2-grad99)/h c print *,ff2 maxstep = 2.0d0; do 200 i = 1,og kkk = maxstep/(2.0d0**i); c print *,'step ',kkk x1aktuell = x1 - kkk*ff1; x2aktuell = x2 - kkk*ff2; if (i == 1) then call bundle(x1aktuell,x2aktuell,ffmin,zustand); x1best = x1aktuell x2best = x2aktuell kkkgenommen = kkk else call bundle(x1aktuell,x2aktuell,ff,zustand); if (ff < ffmin) then ffmin = ff x1best = x1aktuell x2best = x2aktuell kkkgenommen = kkk endif endif 200 continue x1 = x1best x2 = x2best c feld(lala) = kkkgenommen c print *,lala,kkkgenommen if (mod(lala,ausgabeiter) == 0) then call cpu_time(t1) print *,' ' print *,'Iter = ',lala print *,'x1 = ',x1 print *,'x2 = ',x2 print *,'ff = ',ffmin print *,'Schrittweite = ',kkkgenommen print *,'Zeit = ',t1 endif 100 continue c call ausgabe(feld,maxiter) call cpu_time(t1) print *,'Zeit ',t1 end c *********** c c *********** subroutine bundle(x1,x2,ff,zustand) implicit none double precision ff1,ff2,ff3,x1,x2,ff integer zustand ff1 = x1**2.0d0+x2**4.0d0 ff2 = (2.0d0-x1)**2+(2.0d0-x2)**2.0d0 ff3 = 2.0d0*exp(-x1+x2) ff = ff1 zustand = 1 if (ff2 > ff) then ff = ff2 zustand = 2 endif if (ff3 > ff) then ff = ff3 zustand = 3 endif c ff = max(ff1,ff2) c ff = max(ff,ff3) end c ********** c c ********** subroutine subgrad99(x1,x2,ff,zustand,nr) implicit none double precision x1,x2,ff integer zustand,nr if (zustand == 1) then c Ableitung if (nr == 1) then ff = 2.0d0*x1 endif if (nr == 2) then ff = 4.0d0*x2**3.0d0 endif endif if (zustand == 2) then if (nr == 1) then ff = 2.0d0*x1-4.0d0 endif if (nr == 2) then ff = 2.0d0*x2-4.0d0 endif endif if (zustand == 3) then if (nr == 1) then ff = -2.0d0*exp(-x1+x2) endif if (nr == 2) then ff = 2.0d0*exp(-x1+x2) endif endif end c *********** c c *********** subroutine ausgabe(feld,maxiter) implicit none double precision feld(1000000) integer i,maxiter c open(11,file='kurven.dat') do 100 i = 1,maxiter write(11,*) feld(i) 100 continue close(11) return end \sourceoff \showoff


   Profil
hyperG
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 03.02.2017
Mitteilungen: 1923
  Beitrag No.12, eingetragen 2022-10-26

Hallo Delastelle, 130.000 Iterationen und im Ergebnis nur 4 richtige Nachkommastellen...? Das ist doch hochgradig uneffektiv! Um auszuschließen, dass viel Zeit bis zum Finden der 1. richtigen Nachkommastelle "verbraten" wird, folgende Frage: Wie verhalten sich beide von Dir vorgestellten Algorithmen (SQP & Subgrad...), wenn der Startwert bereits mit der 1. Nachkommastelle übereinstimmt (x1=1.1 , x2=0.8)? Da Du nur double verwendest, könnte man dann fragen: Wie viele Iterationen sind vom "guten Startwert" noch nötig, um mindestens 8 richtige Nachkommastellen für beide Algorithmen zu bekommen? Wenn man nun "primitiv linear denkt", könnte man meinen, dass bei einer Iterationszahl von z.B. etwa 1000 eine Genauigkeitssteigerung von (8-1)=7 Stellen pro 1000 Iterationen möglich sein müsste. Dass es hier jedoch nicht linear zugeht, zeigt die Tatsache, dass ich es nicht schaffte mit Mathematica mehr als 14 richtige Nachkommastellen mit numerischen Algorithmen (egal welche Methode) zu bekommen, auch wenn ich über 110 Nachkommastellen genau rechnete! Also untersuchte ich den Einfluss der Exp-Funktion innerhalb dieses Fundbereiches: kein Einfluss, da mit Ergebnis um 1.57 weit unterhalb der 1.95! Also vereinfacht sich die Aufgabe für diesen Spezialfall hier zu Minimize f(a,b)=Max[(2 - a)^2 + (2 - b)^2, a^2 + b^4) im Fundbereich. Dort findet Mathematica nun eine symbolische Lösung für a: \showon Nullstelle einer Funktion 7. Grades, die in eine Nullstellenfunktion 8. Grades hineinverpflanzt wird. Per Newton-Verf. geht das super schnell. Zwar könnte man b auch per Newton berechnen, aber da gibt es bereits eine explizite Lösung für die 2. Variable b {also x2; a & b habe ich nur verwendet, damit die MEGA langen Formeln kürzer werden!}: https://matheplanet.com/matheplanet/nuke/html/uploads/b/47407_MegaWurzel.png \showoff Alle 81 Nachkomma-Stellen des Endergebnisses {also lokales Minimum der gegebenen Max-Funktion} aus Beitrag 9 scheinen also zu stimmen. Auch die Probe mit anderen Punkten dicht um die Fundstelle herum zeigten alle ein schlechteres (größeres) Minimum. Ohne Mathematica hätte ich diese Mega Formeln nicht gefunden. Falls es eine "viel einfachere" Lösung gibt, würde mich das interessieren! Was ich nicht verstehe: Warum kann Mathematica numerisch nicht mehr als 14 Nachkommastellen? Selbst eine primitive https://de.wikipedia.org/wiki/Bisektion sollte doch locker (fast linear wenn der Knick nicht wäre) mehr Stellen schaffen? Liegt es tatsächlich am "Max-Knick", der jede Suche durcheinander bringt? Grüße [Die Antwort wurde nach Beitrag No.10 begonnen.]


   Profil
hyperG
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 03.02.2017
Mitteilungen: 1923
  Beitrag No.13, eingetragen 2022-10-26

Diese Max-Funktion zu minimieren gefällt mir immer besser! In den 3 bekanntesten Suchdatenbanken OEIS, Wolfram, http://wayback.cecm.sfu.ca/ sind alle 3 Konstanten (a, b, Min) NICHT zu finden! Numerische Verfahren beißen sich ab 14 Nachkommastellen die Zähne aus! Richtig schöne Knobelaufgabe für Nachkommastellen oberhalb der 20 Nachkommastellen, wenn man kein hochwertiges https://de.wikipedia.org/wiki/Computeralgebrasystem hat. Wäre auch interessant, was andere wie Derive , Macsyma , Magma , Maple , Mathcad, ... dazu sagen? Könnte ja sein, dass es bereits "Spezialisten" unter den Minimum Findern bei dieser Max-Funktion gibt...


   Profil
Delastelle
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 17.11.2006
Mitteilungen: 2209
  Beitrag No.14, vom Themenstarter, eingetragen 2022-10-26

Hallo hyperG! Schön, dass Dir die Aufgabe gefällt. Interessant in der mathematischen Optimierung sind oft gut funktionierende Beispiele und Verfahren. Manchmal hat man aber auch Beispiele die schlecht zu den Verfahren passen. Du kannst gerne versuchen, mehr Nachkommastellen der Lösung zu finden! In der mathematischen Optimierung hat man oft künstlich gestellte Beispiele (wie hier) oder aber Probleme aus der Praxis. Bei Praxis-Problemen ist oft eine hohe Genauigkeit gar nicht möglich. Wenn ich einen LKW belade kann ich wohl kaum eine höhere Genauigkeit als Gramm oder was eine Waage hergibt nehmen. Für das Beispiel von oben sollten Bundle-Verfahren gut funktionieren. Ich habe aber aktuell keine Implementierung davon vorliegen. Eine Idee: wenn man für Nelder-Mead-Verfahren (ohne Ableitungen) die Genauigkeit und die Iterationszahl vergrößern könnte - würde ich noch mal dieses Verfahren testen. Viele Grüße Ronald


   Profil
hyperG
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 03.02.2017
Mitteilungen: 1923
  Beitrag No.15, eingetragen 2022-10-26

\quoteon(2022-10-26 23:31 - Delastelle in Beitrag No. 14) ... Du kannst gerne versuchen, mehr Nachkommastellen der Lösung zu finden! ... \quoteoff Mit der von mir ( & Mathematica) gefundenen symbolischen Lösung, schaffe ich 1000 Nachkommastellen in nur 0,03 s = 30 ms! Mich würde jedoch interessieren, ob es einen numerischen Algorithmus gibt, der auch nur ansatzweise in diese Genauigkeits- & Zeit-Dimension vorstoßen kann? Deshalb habe ich auch die Fundstelle noch nicht genau angegeben, da man dann ja einen besseren Startwert hätte (und mit Tricks nachträglich optimieren könnte!). Außerdem kann ja auch bei mir (und/oder Mathematica) ein Denkfehler vorliegen... -> ich suche also eine Art Validierung meines Ergebnisses. Grüße


   Profil
Delastelle
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 17.11.2006
Mitteilungen: 2209
  Beitrag No.16, vom Themenstarter, eingetragen 2022-10-26

Hallo, \sourceon Octave % x0 = % 1.1390 % 0.8996 % v = 1.9522 % nerv = 188 x0 = [1; -0.1]; [x0,v,nerv] = nelder_mead_min('bundle',x0) \sourceoff Eine (außer SQP-Verfahren) schnelle Lösung liefert das Nelder-Mead-Verfahren unter Octave. Ich habe mir gerade eine kostenlose optim-Toolbox heruntergeladen. Dort ist ein Nelder-Mead enthalten. ff = 1.952224494 (?) Oben nerv könnten die Iterationen sein. Das wären 188. Edit: nerv sind laut Quelltext die Anzahl der Funktionsauswertungen. Viele Grüße Ronald


   Profil
hyperG
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 03.02.2017
Mitteilungen: 1923
  Beitrag No.17, eingetragen 2022-10-27

\quoteon(2022-10-26 23:31 - Delastelle in Beitrag No. 14) ... Eine Idee: wenn man für Nelder-Mead-Verfahren (ohne Ableitungen) die Genauigkeit und die Iterationszahl vergrößern könnte - würde ich noch mal dieses Verfahren testen. Viele Grüße Ronald \quoteoff Ja, die Idee mit NelderMead ist sehr gut: - starke Einschränkung des Suchbereiches (6 Stellen genauer Startwert) - ohne den exp-Anteil - 100 Nachkommastellen für Zwischenergebnisse - 10^7 Iterationen ergibt schon mal in 1,125 s 51 richtige Nachkommastellen! Die Richtung stimmt! [Die Antwort wurde nach Beitrag No.15 begonnen.]


   Profil
hyperG
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 03.02.2017
Mitteilungen: 1923
  Beitrag No.18, eingetragen 2022-10-27

Validierung des symbolischen Endergebnisses aus Beitrag No.9 mit dem numerischen NelderMead-Algorithmus erfolgreich: - 29 Stellen genauer Startwert/Suchbereich - 190 Nachkommastellen für Zwischenergebnisse - 10^8 Iterationen Obergrenze ergibt schon mal in 1,2 s 96 richtige Nachkommastellen im Ergebnis! ABER a & b {x1 & x2} stimmen nur mit 47 richtigen Nachkommastellen überein! Und genau das macht die Suche bei dieser 2D-Funktion für numerische Alg. so kompliziert: a ist etwas kleiner & b etwas größer als die perfekten Werte und so sucht man in einem Bereich, der etwa die Hälfte der Nachkommastellen des Endergebnisses entspricht! Oder mit anderen Worten: das Tal ist so flach, dass man kaum noch Änderungen im Funktionsergebnis hat, um weitere Suchrichtungen zu bestimmen.


   Profil
hyperG
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 03.02.2017
Mitteilungen: 1923
  Beitrag No.19, eingetragen 2022-10-27

Weitere numerische Ergebnisse und die Abweichung der Ergebnisse zum Idealwert in Nachkommastellen bei einem Suchbereich/Startwert, der schon um die 30 Stellen genau ist: \sourceon nameDerSprache WorkingPrecision | NK f(x,y) | NK x & y 100 | 51 | 25 190 | 96 | 47 290 | 145 | 74 \sourceoff In Worten: um N richtige Nachkommastellen für die Koordinaten {x,y}={a,b}={x1,x2} zu bekommen, muss man bei Mathematicas NelderMead etwa mit 4*N Nachkommastellen genau rechnen, wobei das Funktionsergebnis dabei nur um etwa 2*N Nachkommastellen vom Idealwert abweicht. Es könnte auch sein, dass Mathematica seine Abbruchbedingung grob so auslegt, damit die Rechenzeit nicht zu stark ansteigt {viel Zeit für kaum 2 Nachkommastellen mehr; und auch Gefahren der Endlosschleife}. Und zur Frage in der Überschrift: ja, wenn man einfach grob drauf losrechnet gibt es in der Tat viele numerische Probleme. Besser man fangt erst grob an, um einen besseren Startwert zu bekommen. Dann kann man den Suchbereich immer weiter einschränken, bis die 3. Funktion EXP keinen Einfluss mehr hat. Im letzten Schritt kann man dann entweder numerisch hoch optimieren, oder sogar besser noch pseudo-symbolisch {denn die Polynome 7. und 8. Grades sind ja in Wirklichkeit auch numerische Iterationen - jedoch sehr schnell} zu einer schnellen Lösung kommen, bei der 1000 Nachkommastellen in 30 ms kein Problem sind. Mir gefallen solche Aufgaben, wo man mit 2 völlig unterschiedlichen Algorithmen-Übergruppen {numerisch/pseudo-symbolisch} zum selben Ergebnis kommt... bis in den 1000 Nachkommastellen-Bereich. Grüße Gerd


   Profil
hyperG
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 03.02.2017
Mitteilungen: 1923
  Beitrag No.20, eingetragen 2022-10-28

Da Mathematicas Plot3D beim Hineinzoomen auf 6 Nachkommastellen richtig hässlich wird (auch WorkingPrecision-> 22, PlotPoints -> 1000 & MaxRecursion halfen nicht), habe ich das mal mit meinem https://www.lamprechts.de/gerd/3D-online-Plotter.htm probiert: https://matheplanet.com/matheplanet/nuke/html/uploads/b/47407_MaxPow_3D_Zoom8Stellen.png Beim noch weiter hineinzoomen kommt man an die Grenzen von der double Genauigkeit (x:1.139037...1.139038, y: 0.899559...0.899560) und Grafik-Zoom: https://matheplanet.com/matheplanet/nuke/html/uploads/b/47407_MaxPow_3D_Zoom9Stellen.png Schön zu erkennen das breite flache Band mit kaum unterscheidbarer Höhe. Mit double ist hier ein Minimum nicht genauer bestimmbar.


   Profil
Delastelle
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 17.11.2006
Mitteilungen: 2209
  Beitrag No.21, vom Themenstarter, eingetragen 2022-10-28

Schöne Wolken in Beitrag 20 Bild 2! Ich hatte mal mit Nelder-Mead keine guten Ergebnisse bei vielen Variablen - da waren gradientenähnliche Verfahren deutlich besser. Aber hier mit 2 Variablen ist Nelder-Mead gut. Es kommt natürlich immer sehr auf die Zielfunktion an.


   Profil
Delastelle
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 17.11.2006
Mitteilungen: 2209
  Beitrag No.22, vom Themenstarter, eingetragen 2022-10-30

Ich habe jetzt im Internet noch ein Nelder-Mead in Fortran gefunden. (Es gibt dort auch Nelder-Mead in C und Java.) Es benutzt single precision statt double precision. Eine Zeitmessung habe ich hinzugefügt. Eine Lösung wird nach 20 Iterationen in 0.0468 Sekunden gefunden. The minimum was found at 1.15380 0.886542 The value at the minimum is 1.9558 The number of function evaluations was 41 The number of iterations was 20 Zeit 4.680030E-02 Edit: mit Genauigkeit 1e-6 statt 1e-4 wie oben erhalte ich: The minimum was found at 1.13888 0.899681 The value at the minimum is 1.9522 The number of function evaluations was 81 The number of iterations was 42 Zeit 3.120020E-02 /Edit \showon \sourceon Fortran Program FROSEN c http://www.mikehutt.com/frosen.f c Start -1.9,2 ! Michael F. Hutt ! http://www.mikehutt.com ! Mar. 17, 1998 ! $Id: frosen.f,v 1.4 2007/07/10 12:45:32 mike Exp $ ! ! This program will attempt to minimize Rosenbrock's function using the ! Nelder-Mead simplex method. The program was originally developed in C. ! To be consistent with the way arrays are handled in C, all arrays will ! start from 0. ! ! to compile this program with g77 use: ! g77 -ffree-form -o frosen.exe frosen.f ! * Copyright (c) 1998-2004 ! * ! * Permission is hereby granted, free of charge, to any person obtaining ! * a copy of this software and associated documentation files (the ! * "Software"), to deal in the Software without restriction, including ! * without limitation the rights to use, copy, modify, merge, publish, ! * distribute, sublicense, and/or sell copies of the Software, and to ! * permit persons to whom the Software is furnished to do so, subject to ! * the following conditions: ! * ! * The above copyright notice and this permission notice shall be ! * included in all copies or substantial portions of the Software. ! * ! * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, ! * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF ! * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND ! * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE ! * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION ! * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION ! * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. ! ====================================================================== ! Start of main program Real sp(0:1) Integer d Real e Real scale real t0,t1 call cpu_time(t0) c Rose -1- c sp(0) = -1.2 c sp(1) = 1.0 c Rose -2- c sp(0) = -1.9 c sp(1) = 2.0 c max von 3 Funktionen sp(0) = 1 sp(1) = -0.1 d = 2 e = 1.0e-4 scale = 1.0 CALL SIMPLEX(sp,d,e,scale) call cpu_time(t1) print *,'Zeit ',t1 STOP END ! End of main program ! ====================================================================== ! This is the function to be minimized FUNCTION FUNC(x) c implicit none real x(0:1) double precision ff1,ff2,ff3,ff ff1 = x(0)**2.0d0+x(1)**4.0d0 ff2 = (2.0d0-x(0))**2+(2.0d0-x(1))**2.0d0 ff3 = 2.0d0*exp(-x(0)+x(1)) ff = max(ff1,ff2) ff = max(ff,ff3) FUNC = ff c Rose c FUNC = (100*(x(1)-x(0)*x(0))*(x(1)-x(0)*x(0))+(1.0-x(0))* c c (1.0-x(0))) RETURN END ! ====================================================================== ! This is the simplex routine Subroutine SIMPLEX(start, n, EPSILON, scale) Integer n,row Real start(0:n-1) Real EPSILON Real scale ! Define Constants Parameter(MAX_IT = 1000) Parameter(ALPHA=1.0) Parameter(BETA=0.5) Parameter(GAMMA=2.0) Parameter(MATSIZ=10) ! ====================================================================== ! Variable Definitions ! ! Integer vs = vertex with the smallest value ! Integer vh = vertex with next smallest value ! Integer vg = vertex with largest value ! Integer i,j,m,row ! Integer k = track the number of function evaluations ! Integer itr = track the number of iterations ! Real v = holds vertices of simplex ! Real pn,qn = values used to create initial simplex ! Real f = value of function at each vertex ! Real fr = value of function at reflection point ! Real fe = value of function at expansion point ! Real fc = value of function at contraction point ! Real vr = reflection - coordinates ! Real ve = expansion - coordinates ! Real vc = contraction - coordinates ! Real vm = centroid - coordinates ! Real min ! Real fsum,favg,s,cent ! Real vtmp = temporary array passed to FUNC ! ====================================================================== Integer vs,vh,vg Integer i,j,k,itr Real v(0:MATSIZ,0:MATSIZ-1) Real pn,qn Real f(0:MATSIZ) Real fr,fe,fc Real vr(0:MATSIZ-1) Real ve(0:MATSIZ-1) Real vc(0:MATSIZ-1) Real vm(0:MATSIZ-1) Real vtmp(0:MATSIZ-1) Real mini,fsum,favg,cent,s ! Initialize matrices to 0.0 ! Data v/121 * 0.0/ Data f/11 * 0.0/ Data vr/MATSIZ * 0.0/ Data ve/MATSIZ * 0.0/ Data vc/MATSIZ * 0.0/ Data vm/MATSIZ * 0.0/ Data vtmp/MATSIZ * 0.0/ ! create the initial simplex ! assume one of the vertices is 0.0 pn = scale*(sqrt(n+1.)-1.+n)/(n*sqrt(2.)) qn = scale*(sqrt(n+1.)-1.)/(n*sqrt(2.)) Do 10 i=0,n-1 v(0,i) = start(i) 10 Continue Do 20 i=1,n Do 30 j=0,n-1 If (i-1 .EQ. j) Then v(i,j) = pn + start(j) Else v(i,j) = qn + start(j) EndIf 30 Continue 20 Continue ! find the initial function values Do 40 j=0,n ! put coordinates into single dimension array ! to pass it to FUNC Do 45 m=0,n-1 vtmp(m) = v(j,m) 45 Continue f(j) = FUNC(vtmp) 40 Continue ! Print out the initial simplex ! Print out the initial function values Write(*,*) "Initial Values" Write(*,300) ((v(i,j),j=0,n-1),f(i),i=0,n) k = n+1 ! begin main loop of the minimization Do 50 itr=1,MAX_IT ! find the index of the largest value vg = 0 Do 60 j=0,n If (f(j) .GT. f(vg)) Then vg = j EndIf 60 Continue ! find the index of the smallest value vs = 0 Do 70 j=0,n If (f(j) .LT. f(vs)) Then vs = j EndIf 70 Continue ! find the index of the second largest value vh = vs Do 80 j=0,n If ((f(j) .GT. f(vh)).AND.(f(j).LT.f(vg))) Then vh = j EndIf 80 Continue ! calculate the centroid Do 90 j=0,n-1 cent = 0.0 Do 100 m=0,n If (m .NE. vg) Then cent = cent + v(m,j) EndIf 100 Continue vm(j) = cent/n 90 Continue ! reflect vg to new vertex vr Do 110 j=0,n-1 vr(j) = (1+ALPHA)*vm(j) - ALPHA*v(vg,j) 110 Continue fr = FUNC(vr) k = k+1 If ((fr .LE. f(vh)) .AND. (fr .GT. f(vs))) Then Do 120 j=0,n-1 v(vg,j) = vr(j) 120 Continue f(vg) = fr EndIf ! investigate a step further in this direction If (fr .LE. f(vs)) Then Do 130 j=0,n-1 ve(j) = GAMMA*vr(j) + (1-GAMMA)*vm(j) 130 Continue fe = FUNC(ve) k = k+1 ! by making fe < fr as opposed to fe < f(vs), Rosenbrocks function ! takes 63 iterations as opposed to 64 with e = 1.0e-6 If (fe .LT. fr) Then Do 140 j=0,n-1 v(vg,j) = ve(j) 140 Continue f(vg) = fe Else Do 150 j=0,n-1 v(vg,j) = vr(j) 150 Continue f(vg) = fr EndIf EndIf ! check to see if a contraction is necessary If (fr .GT. f(vh)) Then Do 160 j=0,n-1 vc(j) = BETA*v(vg,j) + (1-BETA)*vm(j) 160 Continue fc = FUNC(vc) k = k+1 If (fc .LT. f(vg)) Then Do 170 j=0,n-1 v(vg,j) = vc(j) 170 Continue f(vg) = fc ! at this point the contraction is not successful, ! we must halve the distance from vs to all the ! vertices of the simplex and then continue. ! 10/31/97 - modified C program to account for ! all vertices, Else Do 180 row=0,n If (row .NE. vs) Then Do 190 j=0,n-1 v(row,j) = v(vs,j)+ c (v(row,j)-v(vs,j))/2.0 190 Continue EndIf 180 Continue Do 185 m=0,n-1 vtmp(m) = v(vg,m) 185 Continue f(vg) = FUNC(vtmp) k = k+1 Do 187 m=0,n-1 vtmp(m) = v(vh,m) 187 Continue f(vh) = FUNC(vtmp) k = k+1 EndIf EndIf ! print out the value at each iteration Write(*,*) "Iteration ",itr Write(*,300) ((v(i,j),j=0,n-1),f(i),i=0,n) ! test for convergence fsum = 0.0 Do 200 j=0,n fsum = fsum + f(j) 200 Continue favg = fsum/(n+1.) s = 0.0 Do 210 j=0,n s = s + ((f(j)-favg)**2.)/n 210 Continue s = sqrt(s) If (s .LT. EPSILON) Then GO TO 220 EndIf 50 Continue ! end main loop of the minimization ! find the index of the smallest value 220 vs = 0 Do 230 j=0,n If (f(j) .LT. f(vs)) Then vs = j EndIf 230 Continue ! print out the minimum Do 240 m=0,n-1 vtmp(m) = v(vs,m) 240 Continue mini = FUNC(vtmp) k = k+1 write(*,*)'The minimum was found at ',v(vs,0),v(vs,1) write(*,250)'The value at the minimum is ',mini write(*,*)'The number of function evaluations was',k write(*,*)'The number of iterations was',itr 250 FORMAT(A29,F7.4) 300 FORMAT(F11.6,F11.6,F11.6) RETURN END ! ====================================================================== \sourceoff \showoff


   Profil
hyperG
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 03.02.2017
Mitteilungen: 1923
  Beitrag No.23, eingetragen 2022-10-30

\quoteon(2022-10-30 12:46 - Delastelle in Beitrag No. 22) Ich habe jetzt im Internet noch ein Nelder-Mead in Fortran gefunden. (Es gibt dort auch Nelder-Mead in C und Java.) Es benutzt single precision statt double precision. Eine Zeitmessung habe ich hinzugefügt. Eine Lösung wird nach 20 Iterationen in 0.0468 Sekunden gefunden. The minimum was found at 1.15380 0.886542 The value at the minimum is 1.9558 The number of function evaluations was 41 The number of iterations was 20 Zeit 4.680030E-02 ... \quoteoff "minimum was found at 1.15380 0.886542" Bedeutet in Wirklichkeit: nur 1 richtige Nachkommastelle gefunden! Das wäre in meinen Augen selbst mit "Rechenschieber" kein "Ergebnis gefunden", sondern grob geschätzt. (da ist ja meine 3D Grafik mit ZOOM 5 Stellen genauer) Den Code in c habe ich gefunden. Dort wird mit normalem double (etwa 15 Nachkommastellen) gerechnet. Vermutlich muss die Abbruchbedingung noch optimiert werden... ( e = 1.0e-4 -> EPSILON ist viel zu groß bei double!)


   Profil
Delastelle hat die Antworten auf ihre/seine Frage gesehen.
Delastelle hat selbst das Ok-Häkchen gesetzt.

Wechsel in ein anderes Forum:
 Suchen    
 
All logos and trademarks in this site are property of their respective owner. The comments are property of their posters, all the rest © 2001-2023 by Matroids Matheplanet
This web site was originally made with PHP-Nuke, a former web portal system written in PHP that seems no longer to be maintained nor supported. PHP-Nuke is Free Software released under the GNU/GPL license.
Ich distanziere mich von rechtswidrigen oder anstößigen Inhalten, die sich trotz aufmerksamer Prüfung hinter hier verwendeten Links verbergen mögen.
Lesen Sie die Nutzungsbedingungen, die Distanzierung, die Datenschutzerklärung und das Impressum.
[Seitenanfang]