83

I use scipy.optimize to minimize a function of 12 arguments.

I started the optimization a while ago and still waiting for results.

Is there a way to force scipy.optimize to display its progress (like how much is already done, what are the current best point)?

Joel Vroom
  • 1,611
  • 1
  • 16
  • 30
Roman
  • 124,451
  • 167
  • 349
  • 456
  • 5
    Have you checked `callback` parameter of your minimization function? – mg007 May 24 '13 at 19:31
  • 1
    For another approach without `callback`, see [Funcgradmon](http://stackoverflow.com/questions/40002172/resuming-an-optimization-in-scipy-optimize/40059852#40059852). It saves all `x f g` values, then can write them to a file for plotting. – denis Oct 15 '16 at 16:35

8 Answers8

52

As mg007 suggested, some of the scipy.optimize routines allow for a callback function (unfortunately leastsq does not permit this at the moment). Below is an example using the "fmin_bfgs" routine where I use a callback function to display the current value of the arguments and the value of the objective function at each iteration.

import numpy as np
from scipy.optimize import fmin_bfgs

Nfeval = 1

def rosen(X): #Rosenbrock function
    return (1.0 - X[0])**2 + 100.0 * (X[1] - X[0]**2)**2 + \
           (1.0 - X[1])**2 + 100.0 * (X[2] - X[1]**2)**2

def callbackF(Xi):
    global Nfeval
    print '{0:4d}   {1: 3.6f}   {2: 3.6f}   {3: 3.6f}   {4: 3.6f}'.format(Nfeval, Xi[0], Xi[1], Xi[2], rosen(Xi))
    Nfeval += 1

print  '{0:4s}   {1:9s}   {2:9s}   {3:9s}   {4:9s}'.format('Iter', ' X1', ' X2', ' X3', 'f(X)')   
x0 = np.array([1.1, 1.1, 1.1], dtype=np.double)
[xopt, fopt, gopt, Bopt, func_calls, grad_calls, warnflg] = \
    fmin_bfgs(rosen, 
              x0, 
              callback=callbackF, 
              maxiter=2000, 
              full_output=True, 
              retall=False)

The output looks like this:

Iter    X1          X2          X3         f(X)      
   1    1.031582    1.062553    1.130971    0.005550
   2    1.031100    1.063194    1.130732    0.004973
   3    1.027805    1.055917    1.114717    0.003927
   4    1.020343    1.040319    1.081299    0.002193
   5    1.005098    1.009236    1.016252    0.000739
   6    1.004867    1.009274    1.017836    0.000197
   7    1.001201    1.002372    1.004708    0.000007
   8    1.000124    1.000249    1.000483    0.000000
   9    0.999999    0.999999    0.999998    0.000000
  10    0.999997    0.999995    0.999989    0.000000
  11    0.999997    0.999995    0.999989    0.000000
Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 11
         Function evaluations: 85
         Gradient evaluations: 17

At least this way you can watch as the optimizer tracks the minimum

Joel Vroom
  • 1,611
  • 1
  • 16
  • 30
  • 21
    this seems super inefficient. you have to call the optimization function again in the callback? does adding the callback this way make the optimization go twice as slow? – abcd Jul 05 '16 at 05:48
  • 2
    I tried this and it seems the callback function is not called at all. The code runs the simulation, but callback does not print anything. Should it pass the values of x in my objective function to callbackF(Xi)? – user3015729 Jan 04 '19 at 16:08
17

Following @joel's example, there is a neat and efficient way to do the similar thing. Following example show how can we get rid of global variables, call_back functions and re-evaluating target function multiple times.

import numpy as np
from scipy.optimize import fmin_bfgs

def rosen(X, info): #Rosenbrock function
    res = (1.0 - X[0])**2 + 100.0 * (X[1] - X[0]**2)**2 + \
           (1.0 - X[1])**2 + 100.0 * (X[2] - X[1]**2)**2


    # display information
    if info['Nfeval']%100 == 0:
        print '{0:4d}   {1: 3.6f}   {2: 3.6f}   {3: 3.6f}   {4: 3.6f}'.format(info['Nfeval'], X[0], X[1], X[2], res)
    info['Nfeval'] += 1
    return res

print  '{0:4s}   {1:9s}   {2:9s}   {3:9s}   {4:9s}'.format('Iter', ' X1', ' X2', ' X3', 'f(X)')   
x0 = np.array([1.1, 1.1, 1.1], dtype=np.double)
[xopt, fopt, gopt, Bopt, func_calls, grad_calls, warnflg] = \
    fmin_bfgs(rosen, 
              x0, 
              args=({'Nfeval':0},), 
              maxiter=1000, 
              full_output=True, 
              retall=False,
              )

This will generate output like 

Iter    X1          X2          X3         f(X)     
   0    1.100000    1.100000    1.100000    2.440000
 100    1.000000    0.999999    0.999998    0.000000
 200    1.000000    0.999999    0.999998    0.000000
 300    1.000000    0.999999    0.999998    0.000000
 400    1.000000    0.999999    0.999998    0.000000
 500    1.000000    0.999999    0.999998    0.000000
Warning: Desired error not necessarily achieved due to precision loss.
         Current function value: 0.000000
         Iterations: 12
         Function evaluations: 502
         Gradient evaluations: 98

However, no free launch, here I used function evaluation times instead of algorithmic iteration times as a counter. Some algorithms may evaluate target function multiple times in a single iteration.

Jinguo Liu
  • 665
  • 6
  • 8
  • While this method is indeed more efficient, the iteration id "Iter" here is inaccurate, because the number of iterations and the number of cost function evaluations are two different things. – liwt31 Dec 17 '22 at 10:19
14

Try using:

options={'disp': True}

to force scipy.optimize.minimize to print intermediate results.

Tom Aranda
  • 5,919
  • 11
  • 35
  • 51
Yan
  • 347
  • 3
  • 3
  • 13
    The documentation suggests that this is the correct answer, but in practice this does not work for me. – bremen_matt Jan 09 '18 at 13:40
  • 24
    The manual says this should be the answer, but as of scipy 1.10 this only outputs information at the end of the minimization -- no the progress of the algorithm or intermediate values. – Juanjo May 31 '18 at 14:09
  • @Juanjo I get your point and you are right it is not printing the progress of the minimization. – muammar Jul 03 '18 at 20:33
  • 1
    did anyone figure out how to get a verbose output? I am also not getting anything after setting `disp: True` in `scipy.optimize.brute`- just the end of the minimization as @Juanjo – tsando Oct 26 '18 at 14:45
  • 2
    This only works at convergence. It is not for printing the intermediate results. – user3015729 Jan 04 '19 at 16:10
  • @user3015729 and others: how are you saving the results? I am getting similar weird behaviour with `differential_evolution`, but it might be related to how I am saving the intermediate results --> https://github.com/scipy/scipy/issues/10325#event-2422335734 – jjrr Jun 19 '19 at 08:55
14

Many of the optimizers in scipy indeed lack verbose output (the 'trust-constr' method of scipy.optimize.minimize being an exception). I faced a similar issue and solved it by creating a wrapper around the objective function and using the callback function. No additional function evaluations are performed here, so this should be an efficient solution.

import numpy as np

class Simulator:
def __init__(self, function):
    self.f = function # actual objective function
    self.num_calls = 0 # how many times f has been called
    self.callback_count = 0 # number of times callback has been called, also measures iteration count
    self.list_calls_inp = [] # input of all calls
    self.list_calls_res = [] # result of all calls
    self.decreasing_list_calls_inp = [] # input of calls that resulted in decrease
    self.decreasing_list_calls_res = [] # result of calls that resulted in decrease
    self.list_callback_inp = [] # only appends inputs on callback, as such they correspond to the iterations
    self.list_callback_res = [] # only appends results on callback, as such they correspond to the iterations

def simulate(self, x, *args):
    """Executes the actual simulation and returns the result, while
    updating the lists too. Pass to optimizer without arguments or
    parentheses."""
    result = self.f(x, *args) # the actual evaluation of the function
    if not self.num_calls: # first call is stored in all lists
        self.decreasing_list_calls_inp.append(x)
        self.decreasing_list_calls_res.append(result)
        self.list_callback_inp.append(x)
        self.list_callback_res.append(result)
    elif result < self.decreasing_list_calls_res[-1]:
        self.decreasing_list_calls_inp.append(x)
        self.decreasing_list_calls_res.append(result)
    self.list_calls_inp.append(x)
    self.list_calls_res.append(result)
    self.num_calls += 1
    return result

def callback(self, xk, *_):
    """Callback function that can be used by optimizers of scipy.optimize.
    The third argument "*_" makes sure that it still works when the
    optimizer calls the callback function with more than one argument. Pass
    to optimizer without arguments or parentheses."""
    s1 = ""
    xk = np.atleast_1d(xk)
    # search backwards in input list for input corresponding to xk
    for i, x in reversed(list(enumerate(self.list_calls_inp))):
        x = np.atleast_1d(x)
        if np.allclose(x, xk):
            break
    
    for comp in xk:
        s1 += f"{comp:10.5e}\t"
    s1 += f"{self.list_calls_res[i]:10.5e}"

    self.list_callback_inp.append(xk)
    self.list_callback_res.append(self.list_calls_res[i])

    if not self.callback_count:
        s0 = ""
        for j, _ in enumerate(xk):
            tmp = f"Comp-{j+1}"
            s0 += f"{tmp:10s}\t"
        s0 += "Objective"
        print(s0)
    print(s1)
    self.callback_count += 1

A simple test can be defined

from scipy.optimize import minimize, rosen
ros_sim = Simulator(rosen)
minimize(ros_sim.simulate, [0, 0], method='BFGS', callback=ros_sim.callback, options={"disp": True})

print(f"Number of calls to Simulator instance {ros_sim.num_calls}")

resulting in:

Comp-1          Comp-2          Objective
1.76348e-01     -1.31390e-07    7.75116e-01
2.85778e-01     4.49433e-02     6.44992e-01
3.14130e-01     9.14198e-02     4.75685e-01
4.26061e-01     1.66413e-01     3.52251e-01
5.47657e-01     2.69948e-01     2.94496e-01
5.59299e-01     3.00400e-01     2.09631e-01
6.49988e-01     4.12880e-01     1.31733e-01
7.29661e-01     5.21348e-01     8.53096e-02
7.97441e-01     6.39950e-01     4.26607e-02
8.43948e-01     7.08872e-01     2.54921e-02
8.73649e-01     7.56823e-01     2.01121e-02
9.05079e-01     8.12892e-01     1.29502e-02
9.38085e-01     8.78276e-01     4.13206e-03
9.73116e-01     9.44072e-01     1.55308e-03
9.86552e-01     9.73498e-01     1.85366e-04
9.99529e-01     9.98598e-01     2.14298e-05
9.99114e-01     9.98178e-01     1.04837e-06
9.99913e-01     9.99825e-01     7.61051e-09
9.99995e-01     9.99989e-01     2.83979e-11
Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 19
         Function evaluations: 96
         Gradient evaluations: 24
Number of calls to Simulator instance 96

Of course this is just a template, it can be adjusted to your needs. It does not provide all information about the status of the optimizer (like e.g. in the Optimization Toolbox of MATLAB), but at least you have some idea of the progress of the optimization.

A similar approach can be found here, without using the callback function. In my approach the callback function is used to print output exactly when the optimizer has finished an iteration, and not every single function call.

Antoine Collet
  • 348
  • 2
  • 14
Henri
  • 348
  • 2
  • 8
  • 3
    Really enjoy your solution. To make it compatible with additionnal `args` for the objective function etc. You can change: `def simulate(self, x, *args)` and `result = self.f(x, *args)` – Antoine Collet Jan 11 '21 at 16:50
3

Which minimization function are you using exactly?

Most of the functions have progress report built, including multiple levels of reports showing exactly the data you want, by using the disp flag (for example see scipy.optimize.fmin_l_bfgs_b).

Bitwise
  • 7,577
  • 6
  • 33
  • 50
2

Below is a solution that works for me :

def f_(x):   # The rosenbrock function
    return (1 - x[0])**2 + 100 * (x[1] - x[0]**2)**2

def conjugate_gradient(x0, f):
    all_x_i = [x0[0]]
    all_y_i = [x0[1]]
    all_f_i = [f(x0)]
    def store(X):
        x, y = X
        all_x_i.append(x)
        all_y_i.append(y)
        all_f_i.append(f(X))
    optimize.minimize(f, x0, method="CG", callback=store, options={"gtol": 1e-12})
    return all_x_i, all_y_i, all_f_i

and by example :

conjugate_gradient([2, -1], f_)

Source

camo
  • 21
  • 4
2

It is also possible to include a simple print() statement in the function to be minimized. If you import the function you can create a wapper.

import numpy as np
from scipy.optimize import minimize


def rosen(X): #Rosenbrock function
    print(X)
    return (1.0 - X[0])**2 + 100.0 * (X[1] - X[0]**2)**2 + \
           (1.0 - X[1])**2 + 100.0 * (X[2] - X[1]**2)**2

x0 = np.array([1.1, 1.1, 1.1], dtype=np.double)
minimize(rosen, 
         x0)
mhass
  • 51
  • 4
  • This would print every function evaluation but not at every iteration. Function evaluation and iterations are different in algorithms such as BFGS. In fact, scipy.optimize.minimize returns the number of iterations and number of function evaluations in two different parameters. – Cristian Arteaga Oct 27 '20 at 21:46
  • This is a solution I considered, but I don't want to print something tens or hundreds of thousands of times. I am not super interested in the exact iteration number so what I thought of as a hack is to only print every time "np.random.randint(1000) == 0" is true. This might add some overhead though. – L. IJspeert Oct 22 '21 at 13:41
0

There you go! (Beware: Most of the time, global variables are bad practice.)

from scipy.optimize import minimize
import numpy as np

f = lambda x, b=.1 : x[0]**2 + b * x[1]**2

x0 = np.array( [2.,2.] )
P = [ x0 ]

def save(x):
    global P
    P.append(x)
    
minimize(f, x0=x0, callback=save)

      fun: 4.608946876190852e-13
 hess_inv: array([[4.99995194e-01, 3.78976566e-04],
       [3.78976566e-04, 4.97011817e+00]])
      jac: array([ 5.42429092e-08, -4.27698767e-07])
  message: 'Optimization terminated successfully.'
     nfev: 24
      nit: 7
     njev: 8
   status: 0
  success: True
        x: array([ 1.96708740e-08, -2.14594442e-06])

print(P)

[array([2., 2.]),
 array([0.99501244, 1.89950125]),
 array([-0.0143533,  1.4353279]),
 array([-0.0511755 ,  1.11283405]),
 array([-0.03556007,  0.39608524]),
 array([-0.00393046, -0.00085631]),
 array([-0.00053407, -0.00042556]),
 array([ 1.96708740e-08, -2.14594442e-06])]
Suuuehgi
  • 4,547
  • 3
  • 27
  • 32