python - 更快的 SciPy 优化

Question

我需要找到具有数千个变量的成本函数的最小值。成本函数只是一个最小二乘法计算，可以使用 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())

Answer 1

如评论中所述，强烈建议为具有 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 倍。最后，我们可以通过设置 jac=True 告诉 minimize 我们的目标函数返回目标和梯度。所以

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

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

PS:您还可以尝试由 cyipopt 接口(interface)的最先进的 Ipopt 求解器包裹。它还提供类似于 scipy.optimize.minimize 的类似 scipy 的界面。