运筹系列21：静态最优化问题解析法

在求解最优化问题$\min f(x)$时，如果使用了目标函数的导数，则称为解析法，否则称之为直接法。直接法包括0.618法、多项式插值法等。本篇介绍解析法，定义
$g(x)=\nabla f(x)$为一阶导。多个y的g(x)按照行进行堆叠后生成的矩阵叫Jacobi矩阵，用J来表示。
$G(x)={\nabla }^{2}f(x)$。构成的矩阵又叫Hessian矩阵。Hessian阵是实对称矩阵，可以特征分解为特征值+正交特征向量。定义Hessian阵的条件数=最大特征值/最小特征值。条件数很大时，会出现震荡问题：

解析法主要分两类：线搜索和置信域搜索两种。

1. 无约束最优化问题

1.1 线搜索方法

线搜索方法基本思路是先确定下降方向d，然后确定步长$\alpha$。时刻记住函数的泰勒展开近似：$f(x+\alpha d)=f(x)+\alpha dg(x)+o(x)$。

• 确定步长的方法包括：
  精确线搜索：规定了方向$d_k$后，在$d_k$上$f(x_k+\alpha d_k)$取得最小值时的$\alpha$称为精确线搜索。求解可以使用黄金分割法、插值法、逆二次插值法、brent法等；
  非精确线搜索：保证每次迭代后$f$下降，但并不保证下降最多。非精确线搜索可以用一阶函数控制变化范围，比如一般要求步长满足Armijo准则作为上界：$f(x_k+\alpha d_k)\le f(x_k)+\rho g_k^T{d}_k\alpha$，此外还要满足如下的准则：Goldstein准则：$f(x_k+\alpha d_k)\ge f(x_k)+(1-\rho )g_k^T{d}_k\alpha$，Wolfe准则：$g(x_k+\alpha d_k{)}^T{d}_k\ge \sigma g_k^T{d}_k$，强Wolfe准则：$|g(x_k+\alpha d_k{)}^T{d}_k|\le -\sigma g_k^T{d}_k$。

• 确定方向的方法包括：
  最速下降法(SD) 的方向$d_k=-g_k/||g_k||$，步长为在方向$d_k$上取最小值。可以理解为用线性方程去拟合曲面，只能确定下降方向，要配合线搜索方法去找到最优点。
  基本Newton法的方向$d_k=-G_k^{-1}{g}_k$，步长为1。可以理解$G_k$度量下的最速下降方向。用二次函数去拟合曲面，可以求得下降最多的点。由于要算二阶导的逆，并且二阶导有可能奇异、不正定，所以一般不用。
  阻尼Newton法的方向$d_k=-G_k^{-1}{g}_k$，步长为在方向$d_k$上取最小值（或使用Wolfe准则求步长）。基本Newton法是拟合曲面的最小值，而阻尼Newton法是优化问题的最小值。
  LM法（Levenberg-Marquardt）用于处理迭代过程中，二阶导不正定、矩阵奇异的情况。LM法的方向$d_k=-(G_k+v_{k}I{)}^{-1}{g}_k$，$v_k$不断乘以2直到$G_k+v_{k}I$正定。LM法得到的方向在Newton方向和负梯度方向之间浮动。
  拟Newton法：用一阶导数计算$H_k$来近似$G_k^{-1}$（或者用$B_k$来近似$G_k$）。拟牛顿法的条件是$H_{k+1}(g_{k+1}-g_k)=x_{k+1}-x_k$，则方向$d_k=-H_k{g}_k$。其中，对称秩1公式是假设$H_{k+1}=H_k+\beta u{u}^T$，称为SR1算法，是一个自对偶的算法，简化表达如下：$\Delta H=\frac{(\Delta x-H\Delta g)(\Delta x-H\Delta g{)}^T}{(\Delta x-H\Delta g)\Delta g}$。对称秩2公式是假设$H_{k+1}=H_k+\beta u{u}^T+\gamma v{v}^T$，称为DFP法，简化表达为：$\Delta H=\frac{\Delta x\Delta x^T}{\Delta x^T\Delta g}-\frac{H\Delta g\Delta g^{T}H}{\Delta g^{T}H\Delta g}$。BFGS公式是针对$B_k$来进行修正的公式，也是秩2算法，是DFP算法的对偶算法，性能更佳。一般取$B_0=I$。DFP和BFGS以一定的比例结合起来的算法称为Broyden族公式。L-BFGS则是不再保存完整的$G_k$近似矩阵，而是保存每一步的变换向量，需要的时候迭代计算出来。
  BB法虽然是负梯度方法，但是步长使用拟牛顿法的条件去求解。
  共轭梯度法使用一阶导信息，迭代产生正交的共轭方向（比如二次优化里面的Gram-Schmidt正交化过程）。共轭梯度的迭代方向为$d_k=-g_k+{\beta }_{k-1}{d}_{k-1}$，其中FR方法的${\beta }_{k-1}=\frac{g_k^T{g}_k}{g_{k-1}^T{g}_{k-1}}$，PRP方法的${\beta }_{k-1}=\frac{g_k^T(g_k-g_{k-1})}{g_{k-1}^T{g}_{k-1}}$，PRP${}^{+}$方法的 β k − 1 = max ⁡ { g k T ( g k − g k − 1 ) g k − 1 T g k − 1 , 0 } \beta_{k-1}=\max\{\frac{g_k^T(g_k-g_{k-1})}{g_{k-1}^Tg_{k-1}},0\} βk−1​=max{gk−1T​gk−1​gkT​(gk​−gk−1​)​,0}。共轭下降法的${\beta }_{k-1}=-\frac{g_k^T{g}_k}{d_{k-1}^T{g}_{k-1}}$，DY法的${\beta }_{k-1}=\frac{g_k^T{g}_k}{d_{k-1}^T(g_k-g_{k-1})}$。

1.2 信赖域方法

信赖域方法是限制在一定的范围内，使用近似函数（例如二阶泰勒展开）同时确定方向和步长；如果合适就加大限制范围，如果不合适就减小限制范围。

1.3 LS问题

接下来是最小二乘问题(LS)求解方法。LS是一类非常重要的问题，下面的方法都是针对这个问题提出的。LS问题主要来源于数据拟合问题：
$\min f(x)=\frac{1}{2}r(x{)}^{T}r(x),x\in R^n$，其中$r(x)$是残差
定义$r(x)$的Jacobi矩阵：$J(x)=[\nabla r_0^T,...,\nabla r_m^T{]}^T$
并令$S(x)={\Sigma }_{i=1}^m{r}_i(x){\nabla }^2{r}_i(x)$，
则$f(x)$的梯度$g(x)=J(x{)}^{T}r(x)$，Hesse阵$G(x)=J(x{)}^{T}J(x)+S(x)$，Newton法求解优化方向方程为$(J_k^T{J}_k+S_k)d_k=-J_k^T{r}_k$。

一般方法是将一阶最优条件转化为一个无约束最优化问题，然后用牛顿法或者置信域法去求解。
Gauss-Newton法用于求解非线性最小二乘（LS）问题。
当$S_k$比较小时，称问题为小剩余量二乘问题，这种情况下忽略$s_k$并使用Newton法迭代的方法称为Gauss-Newton法，即$J_k^T{J}_k{d}_k=-J_k^T{r}_k$。当搜索长度${\alpha }_k=1$时称为基本GN方法；带线搜索的GN方法称为阻尼GN方法；
LS问题也可以使用正交化方法求解。将Newton法方程转化为优化问题$\min ||J_{k}d+r_k|{|}^2$，方法为对$J_k$进行QR分解，求解$R_{k}d=-Q_k^T{r}_k$作为迭代方向。
LMF方法是一种基于LM求解LS问题的方法。类似正交化方法，也是将Newton方程转化为最小化问题，然后使用信赖域方法求解，然后用KKT条件再转化为方程组：$(J_k^T{J}_k+v_{k}I)d_k=-J_k^T{r}_k$，$v_k({\Delta }_k^2-||d_k|{|}^2)=0$，然后看实际函数和近似函数变化的比值，超过一定数值(比如0.75)则扩大范围，低于一定数值(比如0.25)则缩小范围。
上述问题也可以使用Dogleg方法求解，若${\Delta }_k$很小，则使用SD方向与信赖域边界的交点；若${\Delta }_k$很大，则使用GN方向在信赖域内部确定优化点；否则在SD和GN中间选一个方向并求解其与信赖域边界的交点。

3. 约束最优化问题

定义问题：$\min f(x)$
s.t.$g(x)\ge 0$
$h(x)=0$。
一般是用Lagrange+KKT去求解。
Lagrange函数：$L(x,\lambda ,\mu )=f(x)-\lambda g(x)-\mu h(x)$
KKT条件：约束最优化问题最优解的一阶必要条件：${\nabla }_{x}L=0$（定常方程）；$g(x)\ge 0$，$h(x)=0$（原始可行）；$\lambda \ge 0$（对偶可行）；$\lambda g(x)=0$（互补剩余）。最后三个条件可以写作 min ⁡ { λ , g ( x ) } = 0 \min \{\lambda, g(x)\}=0 min{λ,g(x)}=0，求出来的点称为KKT点。对于凸规划问题，KKT点就是全局最优解（凸性相当于二阶最优条件）。

罚函数法也是一种将约束最优化问题转为无约束最优化问题的一种方法。外点罚函数法把约束的平方加到目标函数上，不断提高惩罚因子（比如说 { 1 0 k } \{10^k\} {10k}）直至求解结果进入可行域或者满足一定的终止条件；内点罚函数法又叫障碍函数方法，把约束的倒数（或者约束的负对数）加到目标函数上，不断降低障碍因子直至到达边界或满足一定条件。
乘子罚函数法又叫拉格朗日罚函数法，结合了Lagrange函数的乘子惩罚和外点罚函数法的约束平方惩罚。
序列二次规划（SLSQP）方法是解决一般等式约束最优化问题的方法，Lagrange函数为$L(x,\lambda )=f(x)-\lambda c(x)$。Lagrange-Newton法使用Newton法求解上面函数的KKT稳定点，是一个二次收敛的方法。每一步的迭代问题可以转化为一个二次规划问题（并可以使用拟牛顿法求解），不断迭代进行求解。

4 python代码

4.1 一维搜索

scipy默认使用brent方法，这是结合二分法、割线法和逆二次插值法的一种求根方法。具体原理可以参考这里和这里。

布伦特方法如下：

下面是一个例子：找到$f(x)=\exp [(x-0.7{)}^2]$的最小值。

from scipy import optimize
def f(x):
    return -np.exp(-(x - 0.7)**2)
result = optimize.minimize_scalar(f).x

除了默认的Brent方法，还可以使用Golden（黄金分割法）。

4.2 线搜索方法

首先找到最速下降方向，然后再执行一维搜索。

可以使用
‘Nelder-Mead’，‘Powell’ ，‘CG’ ，‘BFGS’ ，‘Newton-CG’ ，‘L-BFGS-B’，‘TNC’，‘COBYLA’，‘SLSQP’，‘trust-constr’，‘dogleg’，‘trust-ncg’，‘trust-exact’ ，‘trust-krylov’。
下面是个例子：

from scipy.optimize import minimize
def func(x):#目标函数
    result = math.pow(x[0] - x[1], 2) + math.pow(x[0] - 2, 2) + math.pow(x[1] - 3, 4)
    return result
minimize(f, [2, -1])  

默认使用bfgs方法。下面是一些insights:

4.4 约束最优化问题

from scipy.optimize import minimize
func = lambda x:(x[0] - 1)**2 + (x[1] - 2.5)**2
# ineq是≥的意思。
cons = ({'type': 'ineq', 'fun': lambda x:  x[0] - 2 * x[1] + 2},\
        {'type': 'ineq', 'fun': lambda x: -x[0] - 2 * x[1] + 6},\
        {'type': 'ineq', 'fun': lambda x: -x[0] + 2 * x[1] + 2})
bnds = ((0, None), (0, None))
minimize(func, [2, 0], bounds = bnds, constraints=cons)