| |
| Autor |
  Optimierung in 2 Variablen - Numerische Probleme oder nicht? |
|
Delastelle
Senior 
Dabei seit: 17.11.2006
Mitteilungen: 2206
 |  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 
Dabei seit: 03.02.2017
Mitteilungen: 1907
 | 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
Dabei seit: 17.11.2006
Mitteilungen: 2206
| 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
Dabei seit: 03.02.2017
Mitteilungen: 1907
| 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
Dabei seit: 17.11.2006
Mitteilungen: 2206
| 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
Dabei seit: 24.10.2018
Mitteilungen: 4229
 | 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
Dabei seit: 03.02.2017
Mitteilungen: 1907
| 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
Dabei seit: 17.11.2006
Mitteilungen: 2206
| 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
Dabei seit: 03.02.2017
Mitteilungen: 1907
| 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
Dabei seit: 03.02.2017
Mitteilungen: 1907
| 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
Dabei seit: 17.11.2006
Mitteilungen: 2206
| 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
Dabei seit: 17.11.2006
Mitteilungen: 2206
| 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
Dabei seit: 03.02.2017
Mitteilungen: 1907
| 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
Dabei seit: 03.02.2017
Mitteilungen: 1907
| 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
Dabei seit: 17.11.2006
Mitteilungen: 2206
| 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
Dabei seit: 03.02.2017
Mitteilungen: 1907
| 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
Dabei seit: 17.11.2006
Mitteilungen: 2206
| 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
Dabei seit: 03.02.2017
Mitteilungen: 1907
| 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
Dabei seit: 03.02.2017
Mitteilungen: 1907
| 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
Dabei seit: 03.02.2017
Mitteilungen: 1907
| 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
Dabei seit: 03.02.2017
Mitteilungen: 1907
| 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
Dabei seit: 17.11.2006
Mitteilungen: 2206
| 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
Dabei seit: 17.11.2006
Mitteilungen: 2206
| 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
Dabei seit: 03.02.2017
Mitteilungen: 1907
| 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. |
|
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]
|
Home / Seite ohne Frame
Wie unterstütze ich den Matheplaneten mit Geld?
