How to create callable jacobian function for scipy minimize routine?

Question

I can use scipy to perform minimization without using a callable jacobian. I would like to use the callable jacobian, but cannot figure out the correct structure of the output of such a function. I have checked online for examples; the similar questions asked here and here do not have any solutions, and the solution posted here uses importable functions that makes the problem seem unnecessarily complicated. Below is an example that shows a successful minimization without the jacobian, and an unsuccessful attempt at minimizing with the jacobian.

Consider using least squares to find the parameters that provide the best fit of data to an ellipse. First, the data is generated.

import numpy as np
from scipy.optimize import minimize
import matplotlib.pyplot as plt

## generate data
npoints = 1000
prms = (8, 30)
a, b = prms
theta = np.random.uniform(0, 2*np.pi, size=npoints)
_x = a * np.cos(theta)
_y = b * np.sin(theta) 
x = _x + np.random.uniform(-1, 1, size=theta.size) # xpoints with noise
y = _y + np.random.uniform(-1, 1, size=theta.size) # ypoints with noise

## scatter-plot data
fig, ax = plt.subplots()
ax.scatter(x, y, c='b', marker='.')
plt.show()
plt.close(fig)

These are the functions that define the ellipse:

## find parameters of best fit

def f(prms, x, y):
    """ """
    a, b = prms
    return np.square(x/a) + np.square(y/b)

def jacobian(prms, x, y):
    """ """
    a, b = prms
    dfdx = 2 * x / np.square(a)
    dfdy = 2 * y / np.square(b)
    return np.array([dfdx, dfdy])

This is the cost function to minimize:

def error(prms, x, y):
    """ """
    residuals = np.square(1 - f(prms, x, y))
    return np.sum(residuals)

Minimizing without the jacobian:

## try without jacobian
x0 = (10, 20.2)
result = minimize(error, x0=x0, args=(x, y))
# print(result)

The print statement commented out just above outputs the following:

      fun: 11.810701788324323
 hess_inv: array([[ 0.02387587, -0.03327788],
       [-0.03327788,  0.36549839]])
      jac: array([ 4.76837158e-07, -2.38418579e-07])
  message: 'Optimization terminated successfully.'
     nfev: 60
      nit: 12
     njev: 15
   status: 0
  success: True
        x: array([ 8.09464494, 30.03002609])

According to the scipy docs, the function jacobian should return an array of the same size as the input vector. Also, the callable jacobian function is only appropriate for a handful of the methods offered by scipy, including 'L-BFGS-B' and 'SLSQP'. Minimizing with the jacobian:

## try with jacobian 
x0 = (10, 20.2)
result = minimize(error, x0=x0, args=(x, y), method='SLSQP', jac=jacobian)
# print(result)

The following error is raised:

0-th dimension must be fixed to 3 but got 2001

Traceback (most recent call last):
  File "stackoverflow.py", line 39, in <module>
    result = minimize(error, x0=x0, args=(x, y), method='SLSQP', jac=jacobian)
  File "/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/scipy/optimize/_minimize.py", line 611, in minimize
    constraints, callback=callback, **options)
  File "/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/scipy/optimize/slsqp.py", line 426, in _minimize_slsqp
    slsqp(m, meq, x, xl, xu, fx, c, g, a, acc, majiter, mode, w, jw)
_slsqp.error: failed in converting 8th argument `g' of _slsqp.slsqp to C/Fortran array

From the error, I thought the problem was the returned array from calling jacobian contained too many entries, so I tried transposing it. This raised the following error:

0-th dimension must be fixed to 3 but got 2001

Traceback (most recent call last):
  File "stackoverflow.py", line 39, in <module>
    result = minimize(error, x0=x0, args=(x, y), method='SLSQP', jac=jacobian)
  File "/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/scipy/optimize/_minimize.py", line 611, in minimize
    constraints, callback=callback, **options)
  File "/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/scipy/optimize/slsqp.py", line 426, in _minimize_slsqp
    slsqp(m, meq, x, xl, xu, fx, c, g, a, acc, majiter, mode, w, jw)
_slsqp.error: failed in converting 8th argument `g' of _slsqp.slsqp to C/Fortran array

How should the callable function jacobian be output such that it is compatible with the scipy minimize(...) routine?

Answer 1

cannot figure out the correct structure of the output of such a function

do you need to minimize fto a single value? - error use for check_grad. x0 are thetas to your x & y arrays of data!! If this is what you need?

import numpy as np
from scipy.optimize import minimize
import matplotlib.pyplot as plt
from scipy.optimize import check_grad

## generate data
npoints = 1000
prms = (8, 30)
a, b = prms
theta = np.random.uniform(0, 2*np.pi, size=npoints)
_x = a * np.cos(theta)
_y = b * np.sin(theta)
x = _x + np.random.uniform(-1, 1, size=theta.size) # xpoints with noise
y = _y + np.random.uniform(-1, 1, size=theta.size) # ypoints with noise

## scatter-plot data
fig, ax = plt.subplots()
ax.scatter(x, y, c='b', marker='.')
plt.show()
plt.close(fig)

## find parameters of best fit
def f(x: np.ndarray):
    """ """
    a, b = prms
    return np.square(x[0]/a) + np.square(x[1]/b)

def jacobian(x: np.ndarray):
    """ """
    a, b = prms
    dfdx = 2 * x[0] / np.square(a)
    dfdy = 2 * x[1] / np.square(b)
    return np.array([dfdx, dfdy])

x0 = (10, 20.2)
error = check_grad(func=f, grad=jacobian, x0=x0)
assert error < 1e-3

## try with jacobian
x0 = (10, 20.2)
result = minimize(f, x0=x0,  method='SLSQP',   jac=jacobian)
print(result.x)
print(result.fun)
##print(f(result.x))

here your params are being taken from outer scope of functions, not being passed as arguments - I see no problem in such design, as params are constants

P.S. if you need smoothing your circle - use scipy.optimize least_squares or curve_fit or can see LS-circle example