python - 使用 scipy 最小化多变量函数。导数未知

Question

我有一个函数，它实际上是对另一个程序(一些 Fortran 代码)的调用。当我调用此函数 (run_moog) 时，我可以解析 4 个变量，并返回 6 个值。这些值都应该接近于 0(为了最小化)。但是，我将它们组合成这样:np.sum(results**2)。现在我有一个标量函数。我想最小化这个函数，即让 np.sum(results**2) 尽可能接近于零。
注意:当此函数 (run_moog) 接受 4 个输入参数时，它会为依赖于这些参数的 Fortran 代码创建一个输入文件。

我从 the scipy docs 尝试了几种优化方法.但没有一个按预期工作。最小化应该能够限制 4 个变量。这是一个尝试:

from scipy.optimize import minimize # Tried others as well from the docs
x0 = 4435, 3.54, 0.13, 2.4
bounds = [(4000, 6000), (3.00, 4.50), (-0.1, 0.1), (0.0, None)]
a = minimize(fun_mmog, x0, bounds=bounds, method='L-BFGS-B')  # I've tried several different methods here
print a

这给了我

  status: 0
 success: True
    nfev: 5
     fun: 2.3194639999999964
       x: array([  4.43500000e+03,   3.54000000e+00,   1.00000000e-01,
         2.40000000e+00])
 message: 'CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL'
     jac: array([ 0., 0., -54090399.99999981, 0.])
     nit: 0

第三个参数略有变化，其他完全一样。也有 5 次函数调用 (nfev) 但没有迭代 (nit)。 output from scipy is shown here.

Answer 1

几种可能性:

1. 试试 COBYLA。它应该是无衍生的，并支持不平等约束。
2. 您不能通过普通界面使用不同的 epsilon；所以尝试将你的第一个变量缩放 1e4。 (除以进去，乘以回来。)
3. 跳过普通的自动 jacobian 构造函数，并创建您自己的:

假设您正在尝试使用 SLSQP，但您没有提供 jacobian 函数。它为你做了一个。它的代码在 slsqp.py 的 approx_jacobian 中。 .这是一个浓缩版:

def approx_jacobian(x,func,epsilon,*args):
    x0 = asfarray(x)
    f0 = atleast_1d(func(*((x0,)+args)))
    jac = zeros([len(x0),len(f0)])
    dx = zeros(len(x0))
    for i in range(len(x0)):
        dx[i] = epsilon
        jac[i] = (func(*((x0+dx,)+args)) - f0)/epsilon
        dx[i] = 0.0

    return jac.transpose()

您可以尝试将该循环替换为:

    for (i, e) in zip(range(len(x0)), epsilon):
        dx[i] = e
        jac[i] = (func(*((x0+dx,)+args)) - f0)/e
        dx[i] = 0.0

您不能将其作为雅可比矩阵提供给最小化，但为此进行修复很简单:

def construct_jacobian(func,epsilon):
    def jac(x, *args):
        x0 = asfarray(x)
        f0 = atleast_1d(func(*((x0,)+args)))
        jac = zeros([len(x0),len(f0)])
        dx = zeros(len(x0))
        for i in range(len(x0)):
            dx[i] = epsilon
            jac[i] = (func(*((x0+dx,)+args)) - f0)/epsilon
            dx[i] = 0.0

        return jac.transpose()
    return jac

然后您可以像这样调用最小化:

minimize(fun_mmog, x0,
         jac=construct_jacobian(fun_mmog, [1e0, 1e-4, 1e-4, 1e-4]),
         bounds=bounds, method='SLSQP')