如何为 scipy 最小化例程创建可调用的雅可比函数？

Question

我可以使用

scipy

来执行最小化，而无需使用可调用的雅可比矩阵。我想使用可调用的雅可比，但无法弄清楚此类函数输出的正确结构。我在网上查过例子； here和here提出的类似问题没有任何解决方案，而here发布的解决方案使用了可导入的函数，这使得问题看起来不必要地复杂。下面的示例显示了在没有雅可比的情况下成功的最小化，以及在使用雅可比的情况下尝试最小化的失败。

考虑使用最小二乘法来查找为椭圆提供最佳数据拟合的参数。首先，生成数据。

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)

这些是定义椭圆的函数：

## 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])

这是要最小化的成本函数：

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

不使用雅可比最小化：

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

上面注释掉的打印语句输出以下内容：

      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])

根据 scipy 文档，函数

jacobian

应返回与输入向量大小相同的数组。此外，可调用的雅可比函数仅适用于

scipy

提供的少数方法，包括

'L-BFGS-B'

和

'SLSQP'

。使用雅可比最小化：

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

出现以下错误：

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

从错误来看，我认为问题是调用

jacobian

返回的数组包含太多条目，所以我尝试转置它。这引发了以下错误：

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

可调用函数

jacobian

应该如何输出才能与scipy

minimize(...)

例程兼容？

Answer 1

如果这是您所需要的？

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))

您的参数是从函数的外部范围获取的 - 我认为这种设计没有问题，因为

params

是

constants