计算科学 - 如何在 Python 中使用惩罚函数解决优化问题 - 吾爱随笔录

如何在 Python 中使用惩罚函数解决优化问题

计算科学凸优化约束优化 scipy 惩罚法

2021-12-25 19:42:31

我正在实现一个简单的二次优化问题：

min_{x} {\underline{x}}^{T} Q \underline{x}

$\min _x \; {\underline{x}}^T Q {\underline{x}}$

s . t . {\underline{μ}}^{T} \underline{x} = R^{*}

$s.t. \,\quad {\underline{\mu}}^T{\underline{x}} = R^*$

{\underline{1}}^{T} \underline{x} = 1

$\quad \quad \underline{1}^T\underline{x} = 1$ 一旦上述方法起作用，我还希望继续包含不等式约束，作为额外的复杂性。

x_{i} \geq 0

$x_i \geq 0$ 我认为最简单的方法，也是我对实现这些约束最了解的方法，是惩罚函数方法，我们修改目标函数以“引导”优化远离禁区。通过仔细参数化惩罚的大小，我使用 SciPy 的内置 Nelder-Mead Simplex 算法，使用下面的目标函数取得了很好的结果。

def objective(x):
    Q = DF.cov()     # Covariance matrix

    # Penalty Function method

    penalty1 =  0.0005 * abs(np.sum(x)-1)                             # Large for sum(x) <> 1
    penalty2 =  0.05 * abs(R_min - np.matmul(Mus.transpose(), x))     # Large for returns <> R_min

    return np.matmul(x.transpose(),np.matmul(Q,x)) + penalty1 + penalty2

现在，我希望使用其他优化算法（特别是 BFGS 和 Newton-CG），它们需要目标函数的梯度和 Hessian。我已经在无约束情况下实现了导数函数，但是通过将惩罚项添加到目标（以及惩罚的导数到梯度函数）优化失败并出现以下错误：

Warning: Desired error not necessarily achieved due to precision loss.
     Current function value: 0.000056
     Iterations: 0
     Function evaluations: 780
     Gradient evaluations: 96

（以前迭代将是几百）。这严格发生在惩罚 1 中，但不是惩罚 2 本身，所以我的导数对于惩罚 1 是错误的：

penalty1_der = np.sign(x)

或者我不能以这种方式使用L1规范吗？我还尝试用更平滑的二次近似替换约束：

penalty1 = np.matmul((x - vector_ones).transpose(), (x - vector_ones))

但不幸的是，尽管这可以防止错误，但 Minimize() 似乎完全忽略了我的惩罚函数（即使参数大大增加）。

如何实现我的约束，以便我可以使用 BFGS/Newton-CG 解决问题？

1个回答

虽然我同意响应者的普遍共识，即这不是解决问题中最小化问题的最佳方法，但我现在已经解决了挑战，可以回答我自己的问题，以分享在使用惩罚方法时可能克服类似问题的方法解决 Python 中的优化问题。

关键的数学问题确实是惩罚函数的不可微分性；似乎最佳实践是使用与目标函数相同阶的多项式；通过这种方式，您可以确保惩罚的行为与您的目标函数互补。考虑泰勒级数，在概念上似乎很清楚，您应该能够形成这些多项式的总和，具有不同的系数，以很好地逼近您可能想要的几乎任何惩罚。

关键的编程问题（将导致问题中的 SciPy 错误）来自向 scipy.optimize.minimize() 提供不正确的 Jacobian（和/或 Hessian）函数。惩罚函数的加入使梯度向量和 Hessian 矩阵的计算变得更加困难，我不得不手动计算。幸运的是，SciPy 提供了一个函数来测试你的梯度函数：check_grad(F, dF, x_k))，它将你的梯度函数在 x_k 的范数与 x_k 周围小区域的内置有限差分近似进行比较。通常，如果这会返回一些东西 $<10^{-4}$ 那么您的功能可能是正确的（嗯，足够正确）。这不适用于 Hessian 矩阵，因此需要更仔细的计算。

所以，对于这个问题，我用二次近似交换了上面两个惩罚函数中的绝对线性项：

def objective(x):
    Q = DF.cov()                                  # Covariance matrix
    var = np.matmul(W.transpose(),np.matmul(Q,x)) # Variance vector

    # Penalty Function method
    penalty1 =  (np.sum(x)-1)**2                    # Large for sum(x) <> 1
    penalty2 =  100*(R_min - np.dot(Mus, x))**2     # Large for returns <> minR

    return var + penalty1 + penalty2

因此，这些是两次可微的，并且经过一些艰苦的矩阵和向量导数计算，产生了梯度和 Hessian 函数，如下所示：

def der_objective(x):
    Q = DF.cov()

    der = 2 * np.matmul(Q,x)

    penalty1_der = np.array([2*np.sum(x)-2 for i,xi in enumerate(x)])
    penalty2_der = 100*np.array([ 2*Mus[i]*(np.dot(Mus,x) - R_min) for i,xi in enumerate(x)])

    return der + penalty1_der + penalty2_der


def hess_objective(x):
    Q = DF.cov()
    hess = 2 * Q.values

    # Assemble Hessian terms for penalty1 function
    penalty1_hess = 2*np.ones_like(hess)

    # Assemble Hessian terms for penalty2 function
    bcast = np.broadcast(Mus.reshape(len(Mus),1),Mus.reshape(1,len(Mus)))
    penalty2_hess = np.empty(bcast.shape)
    penalty2_hess.flat = 100*np.array([2*a*b for (a,b) in bcast])

    return hess + penalty1_hess + penalty2_hess

然后通常将这些传递给 SciPy 提供的简单接口以选择最小化例程：

result = scipy.optimize.minimize(fun=objective, x0=x_k, method='BFGS',
         jac=der_objective, hess=hess_objective, 
         options={'tol':1e-6,'maxiter':1e3})

其它你可能感兴趣的问题

上一篇固定时间步长的全局截断错误行为下一篇数值分析：切比雪夫系数表示误差