使用scipy最小化多元函数。未知导数。

Question

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

我已经尝试了几种优化方法，来自scipy文档。但没有一种能像预期的那样工作。最小化应该能够在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）。可以在这里查看来自Scipy的输出结果。

- Daniel Thaagaard Andreasen

---

看起来即使没有进行一次迭代，您的函数已经陷入了非常相似的值中。也许像Nelder-Mead这样不使用导数的方法，使用小的ftol和xtol参数可能会解决这个问题。在我看来，网格搜索可能是开始的最佳选择（然后给出初始值以最小化）。您确定不会得到NaN或Inf吗？ （有时即使为参数设置理论界限，算法返回其中之一，如果它不是数值稳定的） - Cristián Antuña

问题在于计算近似梯度的步长太小，导致雅可比矩阵中存在零元素。尝试使用method='L-BFGS-B'和options={'epsilon': 1e-4}，或者一些更大的值（默认为1e-8），直到雅可比矩阵不再有零元素。 - rth

两种解决方案都似乎有所帮助。然而，我想使用其中一个可以使用约束（BFGS、L-BFGS-B、SLSQP）的三个解法之一。所以通过设置 eps: 1e0 来运行，但是在某些时候我会超出我的约束集。从 OP 添加的唯一内容是 minimize 函数中的 , options={'eps': 1e+0})。 - Daniel Thaagaard Andreasen

此外，三个变量的步长相同并不合理。第一个变量的数量级为1e4，而其余变量的数量级为1。epsilon可以是值列表吗？ - Daniel Thaagaard Andreasen

Answer 1

几种可能性：

1. 尝试COBYLA。它应该是无导数的，并支持不等式约束。
2. 你不能通过正常接口使用不同的ε值；所以尝试将第一个变量缩放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

您不能将此作为雅可比矩阵提供给 minimize，但修复它很简单:

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：

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

- Jay Kominek

Answer 2

听起来你的目标函数没有良好的导数行为。输出中的一行 jac: array([ 0., 0., -54090399.99999981, 0.]) 意味着仅更改第三个变量值是显著的。由于对这个变量的导数几乎无限大，所以函数可能存在问题。这也是为什么第三个变量值最终达到了它的最大值。

我建议你至少在参数空间中的几个点上观察导数。使用SciPy的默认步长1e-8和有限差分法计算它们。 这里提供了如何计算导数的示例。

同时，尝试绘制你的目标函数。例如，保持两个参数不变，让另外两个参数变化。如果函数具有多个局部极小值，你就不应该使用基于梯度的方法，如BFGS。

- jasaarim

---

“问题”在于小步长对其他四个参数没有任何影响。然后突然第三个参数发生了微小的变化，导数变得非常大。 目前的解决方案是将所有四个变量带到相同的数量级。这有所帮助。 制作一个图表是个好主意。我甚至没有想到过。谢谢。 - Daniel Thaagaard Andreasen

Answer 3

奈尔德-米德单纯形法（在上面的评论中由Cristián Antuña推荐）被广泛认为是优化（可能表现不佳）函数且没有导数知识的良好选择（参见《C语言数字演算法》第10章）。

您的问题有两个比较具体的方面。第一个是输入的限制，第二个是缩放问题。以下建议解决这些问题，但您可能需要手动迭代几次才能使事情正常。

输入约束

假设您的输入约束形成一个凸区域（如上面的示例所示，但我想将其概括一些），那么您可以编写一个函数

is_in_bounds(p):
    # Return if p is in the bounds

使用此函数，假设算法希望从点from_移动到点to，其中已知from_在该区域内。然后，以下函数将有效地找到在两点之间的直线上可以继续前进的最远点：

from numpy.linalg import norm

def progress_within_bounds(from_, to, eps):
    """
    from_ -- source (in region)
    to -- target point
    eps -- Eucliedan precision along the line
    """

    if norm(from_, to) < eps:
        return from_
    mid = (from_ + to) / 2
    if is_in_bounds(mid):
        return progress_within_bounds(mid, to, eps)
    return progress_within_bounds(from_, mid, eps)

（请注意，这个函数可以为某些区域进行优化，但这几乎不值得费力，因为它甚至不调用您的原始对象函数，而原始对象函数是代价昂贵的。）

（Nelder-Mead 的一个好处是，该函数执行一系列非常直观的步骤。其中一些点显然会使你离开该区域，但很容易做出修改。这里 Nelder Mead 的实现 做了一些标记在形如 ################################################################## 的行之间的修改。）

import copy

'''
    Pure Python/Numpy implementation of the Nelder-Mead algorithm.
    Reference: https://en.wikipedia.org/wiki/Nelder%E2%80%93Mead_method
'''


def nelder_mead(f, x_start, 
        step=0.1, no_improve_thr=10e-6, no_improv_break=10, max_iter=0,
        alpha = 1., gamma = 2., rho = -0.5, sigma = 0.5):
    '''
        @param f (function): function to optimize, must return a scalar score 
            and operate over a numpy array of the same dimensions as x_start
        @param x_start (numpy array): initial position
        @param step (float): look-around radius in initial step
        @no_improv_thr,  no_improv_break (float, int): break after no_improv_break iterations with 
            an improvement lower than no_improv_thr
        @max_iter (int): always break after this number of iterations.
            Set it to 0 to loop indefinitely.
        @alpha, gamma, rho, sigma (floats): parameters of the algorithm 
            (see Wikipedia page for reference)
    '''

    # init
    dim = len(x_start)
    prev_best = f(x_start)
    no_improv = 0
    res = [[x_start, prev_best]]

    for i in range(dim):
        x = copy.copy(x_start)
        x[i] = x[i] + step
        score = f(x)
        res.append([x, score])

    # simplex iter
    iters = 0
    while 1:
        # order
        res.sort(key = lambda x: x[1])
        best = res[0][1]

        # break after max_iter
        if max_iter and iters >= max_iter:
            return res[0]
        iters += 1

        # break after no_improv_break iterations with no improvement
        print '...best so far:', best

        if best < prev_best - no_improve_thr:
            no_improv = 0
            prev_best = best
        else:
            no_improv += 1

        if no_improv >= no_improv_break:
            return res[0]

        # centroid
        x0 = [0.] * dim
        for tup in res[:-1]:
            for i, c in enumerate(tup[0]):
                x0[i] += c / (len(res)-1)

        # reflection
        xr = x0 + alpha*(x0 - res[-1][0])
        ##################################################################
        ##################################################################
        xr = progress_within_bounds(x0, x0 + alpha*(x0 - res[-1][0]), prog_eps)
        ##################################################################
        ##################################################################
        rscore = f(xr)
        if res[0][1] <= rscore < res[-2][1]:
            del res[-1]
            res.append([xr, rscore])
            continue

        # expansion
        if rscore < res[0][1]:
            xe = x0 + gamma*(x0 - res[-1][0])
            ##################################################################
            ##################################################################
            xe = progress_within_bounds(x0, x0 + gamma*(x0 - res[-1][0]), prog_eps)
            ##################################################################
            ################################################################## 
            escore = f(xe)
            if escore < rscore:
                del res[-1]
                res.append([xe, escore])
                continue
            else:
                del res[-1]
                res.append([xr, rscore])
                continue

        # contraction
        xc = x0 + rho*(x0 - res[-1][0])
        ##################################################################
        ##################################################################
        xc = progress_within_bounds(x0, x0 + rho*(x0 - res[-1][0]), prog_eps)
        ##################################################################
        ################################################################## 
        cscore = f(xc)
        if cscore < res[-1][1]:
            del res[-1]
            res.append([xc, cscore])
            continue

        # reduction
        x1 = res[0][0]
        nres = []
        for tup in res:
            redx = x1 + sigma*(tup[0] - x1)
            score = f(redx)
            nres.append([redx, score])
        res = nres

注意：这个实现是GPL的，你可能会觉得可以接受，也可能不行。尽管如此，从任何伪代码修改NM都非常容易，并且你可能想在任何情况下加入模拟退火算法。

扩展

这是一个更棘手的问题，但jasaarim提出了一个有趣的观点。一旦修改后的NM算法找到一个点，你可能想要在固定一些维度的情况下运行matplotlib.contour，以便查看函数的行为。此时，你可能需要重新调整一个或多个维度，并重新运行修改后的NM。

–

- Ami Tavory

Answer 4

有多难得到梯度的解析表达式？如果有了这个，就可以使用有限差分来近似计算海森矩阵与向量的乘积。然后可以使用其他可用的优化例程。

在 SciPy 中可用的各种优化例程中，名为 TNC（具有截断的牛顿共轭梯度）的例程对问题相关的数值非常稳健。

- Hari

---

在这种情况下，很难得到任何解析表达式。 - Daniel Thaagaard Andreasen