更快的 SciPy 优化答案

【问题标题】：Faster SciPy Optimizations更快的 SciPy 优化
【发布时间】：2021-10-01 01:16:45
【问题描述】：

我需要找到具有数千个变量的成本函数的最小值。成本函数只是一个最小二乘计算，可以使用 numpy 向量化轻松快速地计算。尽管如此，优化仍然需要很长的时间。我的猜测是 SciPy 的最小化器而不是我的成本函数中发生了缓慢的运行时间。如何更改 SciPy 的最小化器的参数以加快运行时间？

示例代码：

import numpy as np
from scipy.optimize import minimize

# random data
x = np.random.randn(100, 75)

# initial weights guess
startingWeights = np.ones(shape=(100, 75))

# random y vector 
y = np.random.randn(100)


def costFunction(weights):
   # reshapes flattened weights into 2d matrix
   weights = np.reshape(weights, newshape=(100, 75))

   # weighted row-wise sum
   weighted = np.sum(x * weights, axis=1)

   # squared residuals
   residualsSquared = (y - weighted) ** 2

   return np.sum(residualsSquared)


result = minimize(costFunction, startingWeights.flatten())

【问题讨论】：

你有你的成本函数和你使用的最小化器的示例代码吗？一些最小化器可能会很慢，例如：Nelder Mead，而有些则更快一些，例如 COBYLA。而且由于这些最小化器中的许多经常调用 FORTRAN 实现或类似的实现，因此求解器基本上与它获得的一样快。不过，调整一些参数可能有助于提高速度，例如最大迭代次数。
除了 cmbfast 的评论：您是否提供了准确的客观梯度和粗麻布？如果不是，则两者都由有限差分近似，这反过来导致对梯度的每个评估的目标进行许多评估。因此，提供精确的梯度和 hessian 可以显着加快优化速度。
查看完整的输出，它应该会告诉您调用目标的次数。
添加了示例代码。我使用了基本模型最小化器，没有任何方法规范。请注意，我的方法受此启发：papers.ssrn.com/sol3/papers.cfm?abstract_id=1971363
我的猜测是最小化 7500 个值需要很多函数调用。如果时间随着值的数量呈指数增长，我不会感到惊讶。

标签： python numpy optimization scipy runtime

【解决方案1】：

正如 cmets 中已经指出的，强烈建议为 N = 100*75 = 7500 变量的大问题提供精确的目标梯度。如果没有提供梯度，它将通过有限差分和approx_derivative 函数来近似。但是，有限差分容易出错并且计算量很大，因为每个梯度评估都需要对目标函数进行2*N 评估（没有缓存）。

这可以通过对目标和近似梯度进行计时来轻松说明：

In [7]: %timeit costFunction(startingWeights.flatten())
23.5 µs ± 2.03 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

In [8]: from scipy.optimize._numdiff import approx_derivative

In [9]: %timeit approx_derivative(costFunction, startingWeights.flatten())
633 ms ± 33.3 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

因此，在我的机器上，每次梯度评估都需要超过半秒的时间！评估梯度的更有效方法是算法微分。使用 JAX 库非常简单：

import jax.numpy as jnp
from jax import jit, value_and_grad

def costFunction(weights):
   # reshapes flattened weights into 2d matrix
   weights = jnp.reshape(weights, newshape=(100, 75))

   # weighted row-wise sum
   weighted = jnp.sum(x * weights, axis=1)

   # squared residuals
   residualsSquared = (y - weighted) ** 2

   return jnp.sum(residualsSquared)

# create the derivatives
obj_and_grad = jit(value_and_grad(costFunction))

在这里，value_and_grad 创建了一个评估目标的函数和梯度并返回两者，即obj_value, grad_values = obj_and_grad(x0)。所以让我们计时这个函数：

In [12]: %timeit obj_and_grad(startingWeights.flatten())
132 µs ± 6.62 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

因此，我们现在评估目标和梯度的速度比以前快了近 5000 倍。最后，我们可以告诉minimize，我们的目标函数通过设置jac=True 返回目标和梯度。所以

minimize(obj_and_grad, x0=startingWeights.flatten(), jac=True)

应该会显着加快优化速度。

PS：您还可以尝试使用cyipopt 包接口的最先进的 Ipopt 求解器。它还提供了类似于 scipy.optimize.minimize 的类 scipy 接口。

【讨论】：