13

Ich habe eine Python-Funktion mit 64 Variablen, und ich habe versucht, es mit L-BFGS-B-Methode in der Minimieren-Funktion zu optimieren, jedoch hat diese Methode eine ziemlich starke Abhängigkeit von der ursprünglichen Schätzung und konnte das globale Minimum nicht finden.Wie findet man ein globales Minimum in der Python-Optimierung mit Grenzen?

Aber ich mochte seine Fähigkeit, Grenzen für die Variablen zu setzen. Gibt es einen Weg/Funktion, um das globale Minimum zu finden, während es Grenzen für die Variablen gibt?

+1

Ich vermute, dass http://math.stackexchange.com/ ist ein geeigneterer Ort, um Fragen wie diese zu stellen. – 9000

+0

Könntest du bitte deine Funktion etwas beschreiben - glatt/gradient/hessisch? Wenn Sie es als eine Summe von Quadraten formulieren können, siehe [scipy-optimize-leasesq-with-bound-constraints] (http://stackoverflow.com/questions/9878558/scipy-optimize-leastsq-wit-bound-constraints) . Siehe auch [scicomp.stackexchange.com/search?q=bfgs](http://scicomp.stackexchange.com/search?q=bfgs). – denis

+0

Ich entwerfe 8 Bezier-Kurven im 3D-Raum mit jeweils 6 Kontrollpunkten, die zu minimierende Funktion ist die Merit-Funktion für diese Kurven, die eine lineare Kombination von 4 verschiedenen Parametern (Längen, Krümmungsradius, Nähe, Höhenordnung) ist aus den Kurven abgeleitet. Bisher habe ich versucht scipy.minimize(), bassishopping, aber ich bin immer noch nicht das globale Minimum zu finden – dilyar

Antwort

3

Einige gesunden Menschenverstand Vorschläge für das Debugging und die Visualisierung jedes Optimierer auf Ihre Funktion:

Sind Ihre Zielfunktion und zumutbar Ihre Zwänge?
Wenn die Zielfunktion eine Summe ist, sagen f() + g(), diejenigen separat für alle x in "fx-opt.nptxt" drucken (unten); Wenn f() 99% der Summe ist und g() 1%, untersuchen.

Constraints: wie viele der Komponenten x_i in xfinal an Grenzen stecken, x_i <= lo_i oder >= hi_i?


Wie holprig ist Ihre Funktion im globalen Maßstab?
Run mit mehreren zufälligen Startpunkte und die Ergebnisse speichern zu analysieren/Grundstück:

title = "%s n %d ntermhess %d nsample %d seed %d" % ( # all params! 
    __file__, n, ntermhess, nsample, seed) 
print title 
... 
np.random.seed(seed) # for reproducible runs 
np.set_printoptions(threshold=100, edgeitems=10, linewidth=100, 
     formatter = dict(float = lambda x: "%.3g" % x)) # float arrays %.3g 

lo, hi = bounds.T # vecs of numbers or +- np.inf 
print "lo:", lo 
print "hi:", hi 

fx = [] # accumulate all the final f, x 
for jsample in range(nsample): 
     # x0 uniformly random in box lo .. hi -- 
    x0 = lo + np.random.uniform(size=n) * (hi - lo) 

    x, f, d = fmin_l_bfgs_b(func, x0, approx_grad=1, 
       m=ntermhess, factr=factr, pgtol=pgtol) 
    print "f: %g x: %s x0: %s" % (f, x, x0) 
    fx.append(np.r_[ f, x ]) 

fx = np.array(fx) # nsample rows, 1 + dim cols 
np.savetxt("fx-opt.nptxt", fx, fmt="%8.3g", header=title) # to analyze/plot 

ffinal = fx[:,0] 
xfinal = fx[:,1:] 
print "final f values, sorted:", np.sort(ffinal) 
jbest = ffinal.argmin() 
print "best x:", xfinal[jbest] 

Wenn einige der ffinal Werte einigermaßen gut aussehen, versuchen, mehr zufällige Startpunkte in der Nähe von jenen - , die rein zufällig sicherlich besser als .

Wenn die x s Kurven sind, oder alles real, plotten Sie die besten paar x0 und xfinal.
(A Faustregel Nsample ~ 5 * d oder 10 * d in d Dimensionen zu langsam, zu viele maxiter/maxeval reduzieren, reduzieren ftol -.? Sie nicht ftol für die Exploration 1e-6 brauchen wie dies.)

Wenn Sie reproduzierbare Ergebnisse wünschen, dann müssen Sie ALLE relevanten Parameter in der title und in abgeleiteten Dateien und Plots auflisten, . Sonst wirst du fragen "wo ist diese kommen ??"


Wie holprig ist Ihre Funktion auf Epsilon-Skala ~ 10^-6?
Methoden, die einen Gradienten nähern manchmal ihre letzte Schätzung zurückkehren, aber wenn nicht:

from scipy.optimize import approx_fprime 
for eps in [1e-3, 1e-6]: 
    grad = approx_fprime(x, func, epsilon=eps) 
    print "approx_fprime eps %g: %s" % (eps, grad) 

Wenn jedoch die Gradientenabschätzung schlecht/holprig ist, bevor der Optimierer beenden, Sie das nicht sehen. Dann müssen Sie alle Zwischen [f, x, approx_fprime] speichern, um sie auch zu beobachten; Einfach in Python - fragen Sie, ob das nicht klar ist.

In einigen Problembereichen ist es üblich, von einem angeblichen zu sichern und neu zu starten. Zum Beispiel, wenn Sie auf einer Landstraße verloren sind, zuerst eine Hauptstraße finden, dann von dort neu starten.


Setzen Sie den exect-Aufruf und Version Info auch in Fragen, z.

x, f, d = fmin_l_bfgs_b(func, x0, approx_grad=1, 
      m=ntermhess, factr=factr, pgtol=pgtol, epsilon=eps) 
      # give values for these -- explicit is better than implicit 
versions: numpy 1.8.0 scipy 0.13.2 python 2.7.6 mac 10.8.3 


Zusammenfassung:
erwarten keine Black-Box-Optimierer auf eine Funktion zu arbeiten, die groß angelegte holprig ist, oder epsilon-Skala holprig, oder beides.
Investieren Sie in Testgerüst, und in den Möglichkeiten zu siehe was der Optimierer tut.

1

Vielen Dank an Ihre ausführliche Antwort, aber wie im ziemlich neu in Python, ich habe nicht so recht, wie Sie den Code mein Programm zu implementieren, aber hier war mein Versuch der Optimierung:

x0=np.array((10, 13, f*2.5, 0.08, 10, f*1.5, 0.06, 20, 
      10, 14, f*2.5, 0.08, 10, f*1.75, 0.07, 20, 
      10, 15, f*2.5, 0.08, 10, f*2, 0.08, 20, 
      10, 16, f*2.5, 0.08, 10, f*2.25, 0.09, 20, 
      10, 17, f*2.5, -0.08, 10, f*2.5, -0.06, 20, 
      10, 18, f*2.5, -0.08, 10, f*2.75,-0.07, 20, 
      10, 19, f*2.5, -0.08, 10, f*3, -0.08, 20, 
      10, 20, f*2.5, -0.08, 10, f*3.25,-0.09, 20)) 

# boundary for each variable, each element in this restricts the corresponding element  above 
bnds=((1,12), (1,35), (0,f*6.75), (-0.1, 0.1),(1,35), (0,f*6.75), (-0.1, 0.1),(13, 35), 
    (1,12), (1,35), (0,f*6.75), (-0.1, 0.1),(1,35), (0,f*6.75), (-0.1, 0.1),(13, 35), 
    (1,12), (1,35), (0,f*6.75), (-0.1, 0.1),(1,35), (0,f*6.75), (-0.1, 0.1),(13, 35), 
    (1,12), (1,35), (0,f*6.75), (-0.1, 0.1),(1,35), (0,f*6.75), (-0.1, 0.1),(13, 35), 
    (1,12), (1,35), (0,f*6.75), (-0.1, 0.1),(1,35), (0,f*6.75), (-0.1, 0.1),(13, 35), 
    (1,12), (1,35), (0,f*6.75), (-0.1, 0.1),(1,35), (0,f*6.75), (-0.1, 0.1),(13, 35), 
    (1,12), (1,35), (0,f*6.75), (-0.1, 0.1),(1,35), (0,f*6.75), (-0.1, 0.1),(13, 35), 
    (1,12), (1,35), (0,f*6.75), (-0.1, 0.1),(1,35), (0,f*6.75), (-0.1, 0.1),(13, 35),) 

from scipy.optimize import basinhopping 
from scipy.optimize import minimize 

merit=a*meritoflength + b*meritofROC + c*meritofproximity +d*(distancetoceiling+distancetofloor)+e*heightorder 
minimizer_kwargs = {"method": "L-BFGS-B", "bounds": bnds, "tol":1e0} 
ret = basinhopping(merit_function, x0, minimizer_kwargs=minimizer_kwargs, niter=10, T=0.01) 

zoom = ret['x'] 
res = minimize(merit_function, zoom, method = 'L-BFGS-B', bounds=bnds, tol=1e-5) 
print res 

das Verdienst Funktion kombiniert x0 mit einigen anderen Werten, um 6 Kontrollpunkte für 8 Kurven zu bilden, berechnet dann deren Längen, Krümmungsradien usw. Sie gibt den endgültigen Vorteil als lineare Kombinationen dieser Parameter mit einigen Gewichten zurück.

ich benutzte basinhopping mit niedrigen Genauigkeiten, um die einige Minima zu finden, dann minimize verwendet, um die Präzision des niedrigsten munimum zu erhöhen.

p.s. die Plattform, auf der ich laufe, ist Enthoght canopy 1.3.0, numpy 1.8.0 scipy 0.13.2 mac 10.8.3

+0

Wie viel "springt" Ihre Funktion durch Becken? wäre nützlich für andere, die dies versuchen könnten. Und eine allgemeine Frage: Sind die Kontrollpunkte auf eine große Box, oder jede in einem eigenen 1/6 th eingeschränkt? Wenn eine große Box, , dann ist der Lösungsraum 6! = 720 mal größer als nötig sein muss. – denis

+0

Wie überprüfe ich, wie viel die Funktion hüpft? Wie für 6 Kontrollpunkte, sagen wir P1, P2, P3, P4, P5, P6. In diesen Punkten sind P1, P6 fest, P2y, P2z, P5y, p5z sind festgelegt, daher sind die Variablen P2x, P3x, P3y, P3z, P4x, P4y, P4z, P5x. Jetzt glaube ich, dass P3 und P4 die gleiche große Box teilen, während P2, P5, sich auf einer einzigen x-Achse bewegen – dilyar

+0

Ich würde eine neue Frage vorschlagen, etwas wie "Wie überprüfe ich, wie eine Funktion springt, für scipy.optimize.basehopping ". Aber Sie bekommen vernünftige Ergebnisse, oder? – denis

21

Dies kann mit scipy.optimize.basinhopping erfolgen. Basinhopping ist eine Funktion, die entwickelt wurde, um die globale Minimum einer Zielfunktion zu finden. Es führt wiederholte Minimierungen mit der Funktion scipy.optimize.minimize durch und nimmt nach jeder Minimierung einen zufälligen Schritt im Koordinatenraum. Basinhopping kann Grenzen weiterhin respektieren, indem einer der Minimierer verwendet wird, die Grenzen implementieren (z. B. L-BFGS-B). Hier ist ein Code, der zeigt, wie diese

# an example function with multiple minima 
def f(x): return x.dot(x) + sin(np.linalg.norm(x) * np.pi) 

# the starting point 
x0 = [10., 10.] 

# the bounds 
xmin = [1., 1.] 
xmax = [11., 11.] 

# rewrite the bounds in the way required by L-BFGS-B 
bounds = [(low, high) for low, high in zip(xmin, xmax)] 

# use method L-BFGS-B because the problem is smooth and bounded 
minimizer_kwargs = dict(method="L-BFGS-B", bounds=bounds) 
res = basinhopping(f, x0, minimizer_kwargs=minimizer_kwargs) 
print res 

Der obige Code zu tun, wird für einen einfachen Fall arbeiten, aber man kann immer noch in einem verbotenen Bereich am Ende, wenn basinhopping zufällige Verschiebung Routine, die Sie dort nimmt. das kann zum Glück, indem man einen benutzerdefinierten Schritt außer Kraft gesetzt werden Routine unter dem Schlüsselwort take_step

class RandomDisplacementBounds(object): 
    """random displacement with bounds""" 
    def __init__(self, xmin, xmax, stepsize=0.5): 
     self.xmin = xmin 
     self.xmax = xmax 
     self.stepsize = stepsize 

    def __call__(self, x): 
     """take a random step but ensure the new position is within the bounds""" 
     while True: 
      # this could be done in a much more clever way, but it will work for example purposes 
      xnew = x + np.random.uniform(-self.stepsize, self.stepsize, np.shape(x)) 
      if np.all(xnew < self.xmax) and np.all(xnew > self.xmin): 
       break 
     return xnew 

# define the new step taking routine and pass it to basinhopping 
take_step = RandomDisplacementBounds(xmin, xmax) 
result = basinhopping(f, x0, niter=100, minimizer_kwargs=minimizer_kwargs, 
         take_step=take_step) 
print result 
+0

Eine Variante ohne Schleife: 'return np.clip (x + Zufall (...), xmin, xmax)'. (Mit der gleichen step_size in allen Dimensionen, so vorskalieren, dass das Begrenzungsfeld ungefähr quadratisch ist.) – denis

+0

Hallo. Ich denke du brauchst self.xmin und self.xmax in __call__ statt xmin und xmax? – CPBL