更快的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())

- John Sorensen

---

你有带有成本函数和最小化器的示例代码吗？一些最小化器可能会很慢，例如：Nelder Mead，而有些则稍微快一些，例如COBYLA。由于许多这些最小化器通常调用FORTRAN实现或类似实现，因此求解器基本上是尽可能快的。不过，调整一些参数可能有助于增加速度，例如最大迭代次数。 - cmbfast

除了cmbfast的评论之外：您是否提供了精确的目标梯度和海森矩阵？如果没有，两者都将通过有限差分来近似，这反过来会导致每个梯度评估需要许多目标函数的评估。因此，提供精确的梯度和海森矩阵可以显着加速优化过程。 - joni

添加了示例代码。我使用了基本模型最小化器，没有指定任何方法。请注意，我的方法受到了这篇文章的启发：https://papers.ssrn.com/sol3/papers.cfm?abstract_id=1971363 - John Sorensen

Scipy 不适用于大型问题。 - Erwin Kalvelagen

对于大型/多维数据集，请参考全局优化的统计学习方法，并以批处理方式并行执行优化。 - JeeyCi

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)

应显著加速优化过程。

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

- joni

---

感谢JAX-自动微分。 - JeeyCi

Answer 2

成本函数只是一个最小二乘计算。不要使用最小化，使用最小二乘 - 如果这正是你所需要的 - 例如：

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import least_squares

x= np.array([0.23, 0.66, 0.93, 1.25, 1.75, 2.03, 2.24, 2.57, 2.87, 2.98])
y= np.array([0.25, -0.27, -1.12, -0.45, 0.28, 0.13, -0.27, 0.26, 0.58, 1.03])    
plt.plot(x,y, 'bv')
plt.show()

def model(theta, t):
    return theta[0]* np.exp(t) + theta[1]*np.log(t) + theta[2]*np.sin(t) + theta[3]*np.cos(t)      #  a ∗ e^x   + b ∗ l n x + c ∗ s i n x + d ∗ c o s x

def residual(theta, t, y):       # f-a
     return (model(theta, t) - y).flatten()

theta = np.array([0.5, 0.5, 0.5, 0.5])

res = least_squares(residual, theta,   args=(x, y),  verbose=1)        # jac=jac,
print(res)

x_test = np.linspace(0, 3)
y_test = model(res.x, x_test)
plt.plot(x,y, 'bo', x_test, y_test, 'k-')
plt.show()

• 正如你所看到的，没有地方放你巨大的 x0=startingWeights= shape(100,75)（同样，对于你的最小化方法，你只需要 x0=shape(75) 作为 coeffs_ 的初始猜测——如果你确实需要的话）——因为 x0 只是优化 init_data（100,75)-xs 相对于 init_coeffs_param（或 theta）...

• 优化过程是在每个数据点上（100个ys中的每个数据点）执行的，即使添加了正则化逻辑

• 如果你的问题的形式化是正确的？...如果你更喜欢在优化之前使用距离度量？...或选择全局优化算法来最小化整个空间？...或为最小化方法手动创建Jacobian_fx（作为参数，就像乔尼回答的那样）？...

• 无论如何，我认为，最好不要将你的凸问题重构为 LS，或者只考虑当 线性-LS-问题转换为凸问题时，如果实际上你需要数据 shape(100,75)) 的全局最优解...

• 再次，我怀疑你的数据点包含这么多（75个）！线性独立维度——正则化可以帮助...——要么改变你的最小化目标，要么应用我提到的任何提示

如果可以的话，可以查看LS-Optimization - 其中描述了一些机会（包括Jacobi方法（例如实现），预处理，共轭梯度用于稀疏矩阵[如果您的数据是空间数据，则使用2范数条件数] - 在最小化中选择更好的步进方向进行旋转等）...

附注：当矩阵变大时，需要矩阵分解（When matrices grow up）。

- JeeyCi

---

确实，在minimize vs. least_squares中有一点差异。 - JeeyCi

使用scipy.optimize.least_squares，你甚至可以计算标准差误差。 - JeeyCi

最小二乘问题的正确形式化与简单的最小化问题相比，是OP所期望的。 - JeeyCi

顺便说一句，过度决定问题 Ax = b - JeeyCi