scipy.minimize -- get cost function vs iteration?

The method least_squares does that with parameter verbose=2. However, it is not a general-purpose minimizer, its purpose to to minimize the sum of squares of the given functions. Example:

least_squares(lambda x: [x[0]*x[1]-6, x[0]+x[1]-5], [0, 0], verbose=2)

For other methods, like minimize, there is no such option. Instead of using callback and re-evaluating the cost function, you may want to add some logging to the function itself. For example, here fun appends the computed values to global variable cost_values:

def fun(x):
  c = x[0]**2 - 2*x[0] + x[1]**4
  cost_values.append(c)
  return c

cost_values = []
minimize(fun, [3, 2])
print(cost_values)

In this example there are 4 similar function values for each iteration step, as the minimization algorithm looks around, computing the approximate Jacobian and/or Hessian. So, print(cost_values[::4]) would be the way to get one value of the cost function per step.

But it's not always 4 values per step (depends on dimension and the method used). So it's better to use a callback function to log the costs after each step. The current cost should be stored in a global variable, so it does not have to be recomputed.

def fun(x):
  global current_cost
  current_cost = x[0]**2 - 2*x[0] + x[1]**4
  return current_cost

def log_cost(x):
  cost_values.append(current_cost)  

cost_values = []
minimize(fun, [3, 2], callback=log_cost)
print(cost_values)

This prints

[3.5058199763814986, -0.2358850818406083, -0.56104822688320077, -0.88774448831043995, -0.96018358963745964, -0.98750765702936738, -0.99588975368993771, -0.99867208501468863, -0.99956795994852465, -0.99985981414137615, -0.99995446605426996, -0.99998521591611178, -0.99999519917089297, -0.99999844105574265, -0.99999949379700426, -0.99999983560485239, -0.99999994662329761, -0.99999998266175671]

I figured out a sort of hack using stdlib features, it uses a "deep" redirect of sys.stdout. Note that this does not work with jupyter since IPython hijacks sys.stdout, which removes the .fileno attribute.

It may be possible to patch Jupyter using a tempfile.SpooledTemporaryFile in this way, removing this issue. I don't know.

I believe because this uses OS-level file descriptors, it is also not threadsafe.

import os
import sys
import tempfile

class forcefully_redirect_stdout(object):
    ''' Forces stdout to be redirected, for both python code and C/C++/Fortran
        or other linked libraries.  Useful for scraping values from e.g. the
        disp option for scipy.optimize.minimize.
    '''
    def __init__(self, to=None):
        ''' Creates a new forcefully_redirect_stdout context manager.

        Args:
            to (`None` or `str`): what to redirect to.  If type(to) is None,
                internally uses a tempfile.SpooledTemporaryFile and returns a UTF-8
                string containing the captured output.  If type(to) is str, opens a
                file at that path and pipes output into it, erasing prior contents.

        Returns:
            `str` if type(to) is None, else returns `None`.

        '''

        # initialize where we will redirect to and a file descriptor for python
        # stdout -- sys.stdout is used by python, while os.fd(1) is used by
        # C/C++/Fortran/etc
        self.to = to
        self.fd = sys.stdout.fileno()
        if self.to is None:
            self.to = tempfile.SpooledTemporaryFile(mode='w+b')
        else:
            self.to = open(to, 'w+b')

        self.old_stdout = os.fdopen(os.dup(self.fd), 'w')
        self.captured = ''

    def __enter__(self):
        self._redirect_stdout(to=self.to)
        return self

    def __exit__(self, *args):
        self._redirect_stdout(to=self.old_stdout)
        self.to.seek(0)
        self.captured = self.to.read().decode('utf-8')
        self.to.close()

    def _redirect_stdout(self, to):
        sys.stdout.close() # implicit flush()
        os.dup2(to.fileno(), self.fd) # fd writes to 'to' file
        sys.stdout = os.fdopen(self.fd, 'w') # Python writes to fd

if __name__ == '__main__':
    import re
    from scipy.optimize import minimize

    def foo(x):
        return 1/(x+0.001)**2 + x

    with forcefully_redirect_stdout() as txt:
        result = minimize(foo, [100], method='L-BFGS-B', options={'disp': True})

    print('this appears before `disp` output')
    print('here''s the output from disp:')
    print(txt.captured)
    lines_with_cost_function_values = \
        re.findall(r'At iterate\s*\d\s*f=\s*-*?\d*.\d*D[+-]\d*', txt.captured)

    fortran_values = [s.split()[-1] for s in lines_with_cost_function_values]
    # fortran uses "D" to denote double and "raw" exp notation,
    # fortran value 3.0000000D+02 is equivalent to
    # python value  3.0000000E+02 with double precision
    python_vals = [float(s.replace('D', 'E')) for s in fortran_values]
    print(python_vals)