给出2000 dim变量时,scipy.optimize.minimize('SLSQP')太慢了

Question

我有一个带有约束和上/下界限的非lenear优化问题,所以scipy我必须使用SLSQP.问题显然不凸.我得到了jacobian的目标和约束函数才能正常工作(结果很好/快到300输入向量).所有功能都经过矢量化和调整,以便快速运行.问题是使用1000+输入向量需要很长时间,虽然我可以看到最小化器不是很多地调用我的函数(目标/约束/渐变)并且似乎在内部花费大部分处理时间.我在某地读到了SLSQP的性能是O(n ^ 3).

是否有更好/更快的SLSQP实现或针对python的此类问题的另一种方法？我尝试使用nlopt并以某种方式返回错误的结果,因为我在scipy中使用了完全相同的函数(使用包装器来适应其方法签名).我也没有使用ipopt与pyipopt包,无法使用ipopt二进制文件与python包装器一起工作.

更新:如果它有帮助,我的输入变量基本上是(x,y)元组的向量或表示坐标的2D表面中的点.有了1000分,我最终得到了2000暗淡的输入向量.我想要优化的函数计算彼此之间的点的最佳位置,同时考虑它们的关系和其他约束.所以问题不是稀疏的.

谢谢...

Answer 1

我们对该模型了解不多，但这里有一些注意事项：

1. SLSQP 确实是为小型（密集）、规模良好的模型而设计的。
2. SLSQP 是本地求解器。它将接受非凸问题，但仅提供局部解决方案。
3. 我怀疑 SLSQP 是否存在这种复杂性。无论如何，它并没有过多说明特定问题的性能。
4. IPOPT 是一个大规模、稀疏内点求解器。它将找到本地解决方案。它可以解决非常大的模型。
5. 像 BARON、ANTIGONE 和 COUENNE 这样的全局求解器可以找到全局解决方案（如果您不耐烦并过早停止，则解决方案的质量会受到限制）。这些求解器（大多数时候）比本地求解器慢。我不知道直接的 Python 链接。
6. 如果您有一个良好的起点，那么本地求解器可能正是您所需要的。使用多启动策略可以帮助找到更好的解决方案（尚未证明全局最优，但您可以确信自己没有找到非常糟糕的局部最优）。
7. Python 框架 PYOMO 提供对许多求解器的访问。但是，您将需要重写模型。PYOMO 具有自动微分功能：无需提供梯度。
8. 为了进行测试，您可以尝试在 AMPL 或 GAMS 中转录模型并通过NEOS在线求解。这将允许您尝试许多求解器。AMPL和GAMS都有自动微分功能。

Answer 2

在我看来,scipy.minimze提供了一个直观的优化界面.我发现加速成本函数(最终是渐变函数)可以为你提供很好的加速.

例如,采用N维Rosenbrock函数:

import numpy as np
from scipy.optimize import minimize


def rosenbrock(x, N):
    out = 0.0
    for i in range(N-1):
        out += 100.0 * (x[i+1] - x[i]**2)**2 + (1 - x[i])**2
    return out

# slow optimize
N = 20
x_0 = - np.ones(N)
%timeit minimize(rosenbrock, x_0, args=(N,), method='SLSQP', options={'maxiter': 1e4})
res = minimize(rosenbrock, x_0, args=(N,), method='SLSQP', options={'maxiter': 1e4})
print(res.message)

优化产生的时间

102 ms ± 1.86 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
Optimization terminated successfully.

现在你可以使用numba加速目标函数,并提供一个天真的函数来计算渐变,如下所示:

from numba import jit, float64, int64

@jit(float64(float64[:], int64), nopython=True, parallel=True)
def fast_rosenbrock(x, N):
    out = 0.0
    for i in range(N-1):
        out += 100.0 * (x[i+1] - x[i]**2)**2 + (1 - x[i])**2
    return out


@jit(float64[:](float64[:], int64), nopython=True, parallel=True)
def fast_jac(x, N):
    h = 1e-9
    jac = np.zeros_like(x)
    f_0 = fast_rosenbrock(x, N)
    for i in range(N):
        x_d = np.copy(x)
        x_d[i] += h
        f_d = fast_rosenbrock(x_d, N)
        jac[i] = (f_d - f_0) / h
    return jac

这基本上只是在目标函数中添加一个装饰器,允许并行计算.现在我们可以再次优化时间:

print('with fast jacobian')
%timeit minimize(fast_rosenbrock, x_0, args=(N,), method='SLSQP', options={'maxiter': 1e4}, jac=fast_jac)
print('without fast jacobian')
%timeit minimize(fast_rosenbrock, x_0, args=(N,), method='SLSQP', options={'maxiter': 1e4})
res = minimize(fast_rosenbrock, x_0, args=(N,), method='SLSQP', options={'maxiter': 1e4}, jac=fast_jac)
print(res.message)

尝试两者,有和没有提供快速雅可比功能.这个输出是:

with fast jacobian
9.67 ms ± 488 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
without fast jacobian
27.2 ms ± 2.4 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
Optimization terminated successfully.

对于小小的努力,这是大约10倍的加速.您可以通过此方式实现的改进在很大程度上取决于成本函数的低效率.我有一个成本函数,其中有几个计算,并且能够获得大约10 ^ 2 - 10 ^ 3的加速.

这种方法的优点是它的功能很少,你可以保持scipy和它漂亮的界面.