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?