4

我可以scipy在不使用可调用雅可比的情况下执行最小化。我想使用可调用的 jacobian，但无法弄清楚这种函数输出的正确结构。我在网上查了一些例子；此处和此处提出的类似问题没有任何解决方案，并且此处发布的解决方案使用了可导入的功能，这使问题看起来不必要地复杂。下面是一个示例，显示了没有雅可比的成功最小化，以及使用雅可比最小化的不成功尝试。

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

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 docs，该函数jacobian应返回与输入向量大小相同的数组。此外，可调用的 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以使其与 scipyminimize(...)例程兼容？

4

0 回答 0