Minimierung einer multivariablen Funktion mit scipy. Derivative nicht bekannt

Question

Ich habe eine Funktion, die eigentlich ein Aufruf an ein anderes Programm ist (einige Fortran-Code). Wenn ich diese Funktion (run_moog) aufrufe, kann ich 4 Variablen analysieren, und es gibt 6 Werte zurück. Diese Werte sollten alle nahe bei 0 liegen (um zu minimieren). Allerdings habe ich sie so kombiniert: np.sum(results**2). Jetzt habe ich eine Skalarfunktion. Ich möchte diese Funktion minimieren, d. H. Die np.sum(results**2) so nahe wie möglich an Null bringen.
Hinweis: Wenn diese Funktion (run_moog) die 4 Eingabeparameter verwendet, erstellt sie eine Eingabedatei für den Fortran-Code, der von diesen Parametern abhängt.

Ich habe mehrere Möglichkeiten versucht, dies zu optimieren von the scipy docs. Aber keiner funktioniert wie erwartet. Die Minimierung sollte Grenzen für die 4 Variablen haben können. Hier ist ein Versuch:

from scipy.optimize import minimize # Tried others as well from the docs 
x0 = 4435, 3.54, 0.13, 2.4 
bounds = [(4000, 6000), (3.00, 4.50), (-0.1, 0.1), (0.0, None)] 
a = minimize(fun_mmog, x0, bounds=bounds, method='L-BFGS-B') # I've tried several different methods here 
print a 

Diese dann gibt mir

status: 0 
success: True 
    nfev: 5 
    fun: 2.3194639999999964 
     x: array([ 4.43500000e+03, 3.54000000e+00, 1.00000000e-01, 
     2.40000000e+00]) 
message: 'CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL' 
    jac: array([ 0., 0., -54090399.99999981, 0.]) 
    nit: 0 

Der dritte Parameter ändert sich leicht, während die anderen sind genau die gleichen. Außerdem gab es 5 Funktionsaufrufe (nfev), aber keine Iterationen (nit). Die output from scipy is shown here.

Answer 1

Paar Möglichkeiten:

1. Try COBYLA. Es sollte derivatfrei sein und Ungleichheitsbeschränkungen unterstützen.
2. Sie können verschiedene epsilons nicht über die normale Schnittstelle verwenden; Versuchen Sie also, Ihre erste Variable um 1e4 zu skalieren. (Teilen Sie es in gehen, multiplizieren wieder heraus kommen.)
3. den normalen automatischen jacobian Konstruktor überspringen, und machen Sie Ihre eigenen:

Sagen Sie versuchen SLSQP zu verwenden, und Sie stellen keine jacobian Funktion. Es macht eins für dich. Der Code dafür ist in approx_jacobian in slsqp.py. Hier ist eine verkürzte Version:

def approx_jacobian(x,func,epsilon,*args): 
    x0 = asfarray(x) 
    f0 = atleast_1d(func(*((x0,)+args))) 
    jac = zeros([len(x0),len(f0)]) 
    dx = zeros(len(x0)) 
    for i in range(len(x0)): 
     dx[i] = epsilon 
     jac[i] = (func(*((x0+dx,)+args)) - f0)/epsilon 
     dx[i] = 0.0 

    return jac.transpose() 

Sie könnten versuchen, diese Schleife zu ersetzen mit:

for (i, e) in zip(range(len(x0)), epsilon): 
     dx[i] = e 
     jac[i] = (func(*((x0+dx,)+args)) - f0)/e 
     dx[i] = 0.0 

Sie können das nicht bieten als jacobian zu minimize, aber es oben für die Fixierung ist einfach:

def construct_jacobian(func,epsilon): 
    def jac(x, *args): 
     x0 = asfarray(x) 
     f0 = atleast_1d(func(*((x0,)+args))) 
     jac = zeros([len(x0),len(f0)]) 
     dx = zeros(len(x0)) 
     for i in range(len(x0)): 
      dx[i] = epsilon 
      jac[i] = (func(*((x0+dx,)+args)) - f0)/epsilon 
      dx[i] = 0.0 

     return jac.transpose() 
    return jac 

können Sie rufen dann minimize wie:

minimize(fun_mmog, x0, 
     jac=construct_jacobian(fun_mmog, [1e0, 1e-4, 1e-4, 1e-4]), 
     bounds=bounds, method='SLSQP') 

Answer 2

Es hört sich an, als ob Ihre Zielfunktion nicht gut behaltende Ableitungen hat. Die Zeile in der Ausgabe jac: array([ 0., 0., -54090399.99999981, 0.]) bedeutet, dass nur der dritte Variablenwert geändert werden kann. Und weil das Derivat w.r.t. zu dieser Variable ist praktisch unendlich, da ist wahrscheinlich etwas falsch in der Funktion. Deshalb endet auch der dritte Variablenwert in seinem Maximum.

Ich würde vorschlagen, dass Sie die Derivate betrachten, zumindest in ein paar Punkten in Ihrem Parameterraum. Berechne sie mit finiten Differenzen und der Standardschrittweite von SciPys fmin_l_bfgs_b, 1e-8. Here ist ein Beispiel dafür, wie Sie die Derivate berechnen könnten.

Versuchen Sie auch, Ihre Zielfunktion zu plotten. Zum Beispiel, halten Sie zwei der Parameter konstant und lassen Sie die beiden anderen variieren. Wenn die Funktion mehrere lokale Optima hat, sollten Sie keine gradientenbasierten Methoden wie BFGS verwenden.

Answer 3

Wie schwierig ist es, einen analytischen Ausdruck für den Gradienten zu erhalten? Wenn Sie das haben, können Sie das Produkt des Hessischen mit einem Vektor unter Verwendung der endlichen Differenz approximieren. Dann können Sie andere verfügbare Optimierungsroutinen verwenden.

Unter den verschiedenen Optimierungsroutinen, die in SciPy zur Verfügung stehen, ist TNC (Newton Conjugate Gradient with Truncation) ziemlich robust gegenüber den numerischen Werten, die mit dem Problem verbunden sind.

Answer 4

Die Nelder-Mead Simplex Method (vorgeschlagen von Cristián Antuña in den Kommentaren oben) bekannt ist auch eine gute Wahl zu sein für die Optimierung (posibly ungezogene) Funktionen ohne Kenntnis von Derivaten (siehe Numerical Recipies In C, Chapter 10).

Es gibt zwei etwas spezifische Aspekte zu Ihrer Frage. Die erste ist die Beschränkung der Eingaben und die zweite ist ein Skalierungsproblem. Im Folgenden werden Lösungen für diese Punkte vorgeschlagen, aber Sie müssen möglicherweise manuell einige Male zwischen ihnen durchlaufen, bis die Dinge funktionieren.

Eingangs Constraints

Ihre Eingabe Zwänge Unter der Annahme, bilden eine convex region (als Beispiele oben zeigen, aber ich möchte es ein wenig verallgemeinern), dann können Sie eine Funktion

is_in_bounds(p): 
    # Return if p is in the bounds 

schreiben

Mit dieser Funktion wird angenommen, dass der Algorithmus vom Punkt from_ zum Punkt to wechseln möchte, von dem bekannt ist, dass from_ in der Region liegt. Dann wird die folgende Funktion effizient den äußersten Punkt auf der Linie finden zwischen den beiden Punkten, auf dem sie gehen können:

from numpy.linalg import norm 

def progress_within_bounds(from_, to, eps): 
    """ 
    from_ -- source (in region) 
    to -- target point 
    eps -- Eucliedan precision along the line 
    """ 

    if norm(from_, to) < eps: 
     return from_ 
    mid = (from_ + to)/2 
    if is_in_bounds(mid): 
     return progress_within_bounds(mid, to, eps) 
    return progress_within_bounds(from_, mid, eps) 

(Beachten Sie, dass diese Funktion für einige Regionen optimiert werden kann, aber es ist kaum lohnt sich die Mühe, als es nennt nicht einmal Ihre ursprüngliche Objektfunktion, die die teure ist.)

Einer der netten Aspekte von Nelder-Mead ist, dass die Funktion eine Reihe von Schritten ausführt, die so intuitiv sind. Einige dieser Punkte können Sie offensichtlich aus der Region werfen, aber es ist einfach, dies zu ändern. Hier ist ein implementation of Nelder Mead mit Modifikationen zwischen Paaren von Zeilen der Form markiert ##################################################################:

import copy 

''' 
    Pure Python/Numpy implementation of the Nelder-Mead algorithm. 
    Reference: https://en.wikipedia.org/wiki/Nelder%E2%80%93Mead_method 
''' 


def nelder_mead(f, x_start, 
     step=0.1, no_improve_thr=10e-6, no_improv_break=10, max_iter=0, 
     alpha = 1., gamma = 2., rho = -0.5, sigma = 0.5): 
    ''' 
     @param f (function): function to optimize, must return a scalar score 
      and operate over a numpy array of the same dimensions as x_start 
     @param x_start (numpy array): initial position 
     @param step (float): look-around radius in initial step 
     @no_improv_thr, no_improv_break (float, int): break after no_improv_break iterations with 
      an improvement lower than no_improv_thr 
     @max_iter (int): always break after this number of iterations. 
      Set it to 0 to loop indefinitely. 
     @alpha, gamma, rho, sigma (floats): parameters of the algorithm 
      (see Wikipedia page for reference) 
    ''' 

    # init 
    dim = len(x_start) 
    prev_best = f(x_start) 
    no_improv = 0 
    res = [[x_start, prev_best]] 

    for i in range(dim): 
     x = copy.copy(x_start) 
     x[i] = x[i] + step 
     score = f(x) 
     res.append([x, score]) 

    # simplex iter 
    iters = 0 
    while 1: 
     # order 
     res.sort(key = lambda x: x[1]) 
     best = res[0][1] 

     # break after max_iter 
     if max_iter and iters >= max_iter: 
      return res[0] 
     iters += 1 

     # break after no_improv_break iterations with no improvement 
     print '...best so far:', best 

     if best < prev_best - no_improve_thr: 
      no_improv = 0 
      prev_best = best 
     else: 
      no_improv += 1 

     if no_improv >= no_improv_break: 
      return res[0] 

     # centroid 
     x0 = [0.] * dim 
     for tup in res[:-1]: 
      for i, c in enumerate(tup[0]): 
       x0[i] += c/(len(res)-1) 

     # reflection 
     xr = x0 + alpha*(x0 - res[-1][0]) 
     ################################################################## 
     ################################################################## 
     xr = progress_within_bounds(x0, x0 + alpha*(x0 - res[-1][0]), prog_eps) 
     ################################################################## 
     ################################################################## 
     rscore = f(xr) 
     if res[0][1] <= rscore < res[-2][1]: 
      del res[-1] 
      res.append([xr, rscore]) 
      continue 

     # expansion 
     if rscore < res[0][1]: 
      xe = x0 + gamma*(x0 - res[-1][0]) 
      ################################################################## 
      ################################################################## 
      xe = progress_within_bounds(x0, x0 + gamma*(x0 - res[-1][0]), prog_eps) 
      ################################################################## 
      ################################################################## 
      escore = f(xe) 
      if escore < rscore: 
       del res[-1] 
       res.append([xe, escore]) 
       continue 
      else: 
       del res[-1] 
       res.append([xr, rscore]) 
       continue 

     # contraction 
     xc = x0 + rho*(x0 - res[-1][0]) 
     ################################################################## 
     ################################################################## 
     xc = progress_within_bounds(x0, x0 + rho*(x0 - res[-1][0]), prog_eps) 
     ################################################################## 
     ################################################################## 
     cscore = f(xc) 
     if cscore < res[-1][1]: 
      del res[-1] 
      res.append([xc, cscore]) 
      continue 

     # reduction 
     x1 = res[0][0] 
     nres = [] 
     for tup in res: 
      redx = x1 + sigma*(tup[0] - x1) 
      score = f(redx) 
      nres.append([redx, score]) 
     res = nres 

Hinweis Diese Implementierung GPL ist, die entweder gut für Sie ist oder nicht. Es ist extrem einfach, NM von irgendeinem Pseudocode zu modifizieren, und Sie könnten simulated annealing auf jeden Fall einwerfen.

Skalierung

Dies ist ein heikler Problem, aber jasaarim hat einen interessanten Punkt in Bezug auf das macht. Sobald der modifizierte NM-Algorithmus einen Punkt gefunden hat, möchten Sie möglicherweise matplotlib.contour ausführen, während Sie einige Dimensionen korrigieren, um zu sehen, wie sich die Funktion verhält. An dieser Stelle möchten Sie möglicherweise eine oder mehrere Dimensionen neu skalieren und die geänderte NM erneut ausführen.

-