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

Question

我有一个带有约束和上限/下限的非线性优化问题，所以对于 scipy，我必须使用 SLSQP。这个问题显然不是凸的。 我让目标函数和约束函数的雅可比函数都能正常工作（结果很好/很快，最多 300 个输入向量）。所有功能都经过矢量化和调整以非常快速地运行。问题是使用 1000+ 输入向量需要很长时间，尽管我可以看到最小化器并没有大量调用我的函数（目标/约束/梯度）并且似乎在内部花费了大部分处理时间。我在某处读到 SLSQP 的性能是 O(n^3)。

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

更新：如果有帮助，我的输入变量基本上是（x，y）元组或表示坐标的二维表面中的点的向量。有了 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

令人惊讶的是，在我更改成本函数以包含不等式约束和边界之后，我发现了一个相对不错的解决方案，它使用深度学习框架 Tensorflow 的优化器，使用基本梯度下降（实际上是 RMSProp，带动量的梯度下降）约束作为惩罚（我想这与拉格朗日方法相同）。它训练速度超快，并在约束惩罚上使用适当的 lambda 参数快速收敛。我什至不必重写 jacobians，因为 TF 会处理这一点，显然不会对速度产生太大影响。

在此之前，我设法让 NLOPT 工作，它比 scipy/SLSQP 快得多，但在更高维度上仍然很慢。此外，NLOPT/AUGLANG 速度非常快，但收敛性很差。

这就是说，在 20k 变量时，它仍然很慢。部分是由于内存交换和成本函数与成对欧几里得距离至少为 O(n^2)（我使用 (x-x.t)^2+(y-y.t)^2 进行广播）。所以仍然不是最优的。

Answer 3

在我看来 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 及其漂亮的界面。