lecture04 约束最优化问题

Lecture 4 约束最优化问题

目录

• 等式约束
  • 一阶必要条件与牛顿法
  • 高斯牛顿法
• 不等式约束
  • 一阶必要条件（KKT条件）
  • Active-Set法
  • 障碍函数法/内点法
  • 罚函数法
  • 增广拉格朗日法
• 二次规划QP

等式约束

考虑如下带等式约束的最优化问题
$$\begin{array}{rl}\underset{x}{\min }f(x)& ;f(x):{\mathbb{R}}^n⟶\mathbb{R}\\ s.t.c(x)=0& ;c(x):{\mathbb{R}}^n⟶{\mathbb{R}}^m\end{array}\text{\hspace{1em}(1)}$$
类似于无约束最优化问题，极值点处的一阶必要条件有：

• 函数$f(x)$在无约束/自由的方向的梯度要为0。$x$不自由的方向则是$c(x)$的梯度方向，因为函数受限于$c(x)=0$，$c(x)$的值是不能发生变化的（梯度方向则是$c(x)$的值变化的方向）
• $c(x)=0$

第一点换句话说，$\nabla f(x)$的方向必须要是不自由的方向，在自由方向不能有分量，因此，$\nabla f$的方向要垂直于约束流形（constraint manifold）的方向。
$$\begin{array}{r}\nabla f+\lambda \nabla c=0\text{for some}\lambda \in \mathbb{R}\end{array}\text{\hspace{1em}(2)}$$
在$x$是向量的情况有，
$$\begin{array}{r}\frac{\partial f}{\partial x}+{\lambda }^T\frac{\partial c}{\partial x}=0\end{array}\text{\hspace{1em}(3)}$$
基于此我们定义拉格朗日函数，包含$x\lambda$两个变量。

$$\begin{array}{r}L(x,\lambda )=f(x)+{\lambda }^{T}c(x)\end{array}\text{\hspace{1em}(4)}$$
现在对这个函数的两个变量求导，正好可以得到
$$\begin{array}{rl}{\nabla }_{x}L(x,\lambda )& =\nabla f+(\frac{\partial c}{\partial x}{)}^T\lambda =0\\ {\nabla }_{\lambda }L(x,\lambda )& =c(x)\end{array}\text{\hspace{1em}(5)}$$
对应上面提到的两个一阶必要条件，这也是等式约束下的KKT条件。

ps：一般的课可能讲到这里就停了，因为对于我在学校学的那些简单函数说，这两个条件列出的所有方程有时候就可以全部解析的求出所有的$x$和$\lambda$了，然后再判断各个点是不是极小值点，得到全局的极小值。但是对于实际许多复杂的非线性函数而言，式(5)无法解析求得，需要用迭代求根的方式来求解局部极小值。

现在要解(5)的方程，使用牛顿法，因为这里有两个变量和两个输出（并不一定二维），对(5)式进行一阶泰勒展开，并令其为0（因为是求根），得到(5)的根的迭代式子，也就是(4)的局部极值点的迭代式子（回忆一下，在Lecture 4中介绍单变量的牛顿法，我们也是那么做的）
$$\begin{array}{rl}{\nabla }_{x}L(x+\Delta x,\lambda +\Delta \lambda )& \approx {\nabla }_{x}L(x,\lambda )+\frac{{\partial }^{2}L}{\partial x^2}\Delta x+\frac{{\partial }^{2}L}{\partial x\partial \lambda }\Delta \lambda \\ {\nabla }_{\lambda }L(x+\Delta x,\lambda )& \approx c(x)+\frac{\partial c}{\partial x}\Delta x\end{array}\text{\hspace{1em}(6)}$$
注意，从(5)的第二条式子中，可以得$\frac{{\partial }^{2}L}{\partial x\partial \lambda }=(\frac{\partial c}{\partial x}{)}^T$,最终，得到$x\lambda$的迭代式
$$\begin{array}{r}\left[\begin{array}{cc}\frac{{\partial }^{2}L}{\partial x^2}& (\frac{\partial c}{\partial x}{)}^T\\ \frac{\partial c}{\partial x}& 0\end{array}\right]\left[\begin{array}{c}\Delta x\\ \Delta \lambda \end{array}\right]=\left[\begin{array}{c}-{\nabla }_{x}L(x,\lambda )\\ -c(x)\end{array}\right]\end{array}\text{\hspace{1em}(7)}$$
前面这个大矩阵正是拉格朗日函数的Hessian矩阵，这个系统称为KKT系统。

ps：因为这个矩阵长得比较特殊也比较通用，因此，教授课上说有许多针对性的求解器来迭代这个方程。

高斯牛顿法

式子(7)中有一项对$L$求二阶导数，而后面那项要对约束求二阶导，在实际的问题中，$f$是期望的目标函数，可以设计地比较规范，而$c$是约束，来自实际的物理系统，往往比较复杂也比较难以微分，因此，在对式子(7)迭代时候，可以考虑丢弃(8)中的第二项，这种方法叫做高斯牛顿法。
$$\begin{array}{r}\frac{{\partial }^{2}L}{\partial x^2}={\nabla }^{2}f+\frac{\partial }{\partial x}[(\frac{\partial c}{\partial x}{)}^T\lambda ]\end{array}\text{\hspace{1em}(8)}$$
高斯牛顿法等价于先对系统线性化后在求解线性化后的系统的极值点。

案例分析

使用牛顿法和高斯牛顿法分别对下图所示的问题进行最小化。其中，一圈一圈的是目标函数$f$的等值线，黄色的抛物线是约束。可以看出，牛顿法在不合理的起始点（-3，2）中，反而有可能无法收敛于极小值点。而高斯牛顿法不存在这个问题，因为高斯牛顿法去掉第二项之后，保证Hessian矩阵是正定的，因此，每一步迭代都是向着下降的方向迭代。虽然从理论上牛顿法在靠近极值点附近收敛更快，但是实际来看，高斯牛顿法表现地更加稳定。

牛顿法 初始点(-1,-1)

牛顿法 初始点（-3，2）

高斯牛顿法 初始点(-1,-1)

高斯牛顿法 初始点（-3，2）

ps：从迭代的过程中来看，在迭代过程中$c(x)$并不是一直被满足的，但是，假如迭代收敛，那么最终一定会在约束上。

小结

• (7)中KKT系统的牛顿法迭代也可能收敛于鞍点或极大值点（即使$f$本身是凸的，但是综合的KKT系统还是可能非凸），因此也需要采用regularization。但是，这种regularization与Lecture 3中提到的无约束regularization还有一些区别，后面讲。
• 高斯牛顿法在实际中往往比较常用，因为每次迭代比较快，而且具有超线性的收敛性。

不等式约束

考虑如下带约束优化问题，
$$\begin{array}{r}\underset{x}{\min }f(x)\\ \ni c(x)\ge 0\end{array}\text{\hspace{1em}(9)}$$
极值点处的一阶必要条件同样是：

• 需要$\nabla f(x)=0$在无约束/自由的方向。
• 满足$c(x)\ge 0$

但是，要从数学上表示这两个条件，不等式的情况比等式的情况要复杂的多。
$$\begin{array}{rl}\nabla f-(\frac{\partial c}{\partial x}{)}^T\lambda & =0\leftarrow \text{"Stationarity"}\\ c(x)& \ge 0\leftarrow \text{"Primal Feasibility"}\\ \lambda & \ge 0\leftarrow \text{"Dual Feasibility"}\\ {\lambda }^{T}c(x)& =0\leftarrow \text{"Complementarity"}\end{array}\text{\hspace{1em}(10)}$$
第四个条件互补松弛条件，表明$\lambda$和$c(x)$至少要一个为0。

• 如果$c(x)$不为0，说明约束$c(x)$不起作用（$x$不在边界上），那么此时$\lambda$必须为0（对应于不为0的那个$c(x)$分量），此时第一个条件则变成无约束问题的极值必要条件。
• 如果$c(x)=0$，那么$\lambda \ge 0$，则条件1与普通等式约束的问题的差别在于，$\nabla f$不为0的分量必须是朝着被约束的那边的（不可行的那边）

从(10)中可以看出，不等式约束的一阶必要条件在数学上是比较复杂的，而且还因为互补松弛条件引入了非常强大的非线性，甚至是一种if-else的逻辑，因此要求解(10)得到极小值点是比较复杂的。我们只能另辟蹊径。

并且，KKT条件只是告诉我们，在极值点处会是一个怎样的情况，并没有说如何迭代或者如何优化才能到极值点，因此具体如何解这个问题，以及各个case怎么处理，要看具体的求解器。

不等式约束的求解方法

这里老师并没有详细讲，抽空详细补充一下。

有效集法（Active set method）

核心想法：找一堆guess heuristic function来不断猜测哪些不等式约束是生效的（等式约束），找到之后，不生效的约束则不作任何处理（当作没有），则按照只有等式约束的问题来迭代一次，之后继续判断哪些集合是生效的。

• 当有一个很好的heuristic的时候，这个是很快的。教授在课上举例说，MIT的某个组在DARPA挑战赛中设计了一堆合理的启发函数来估计active set，进而解机器人QP问题，这个启发函数在95%的情况都是正确的，在剩下的5%的情况就再继续估计，这样做可以做到几千赫兹
• 当估计的很差的时候，这种解法效果不好

障碍函数法/内点法（Barrier / Interior-Point Method）

**核心想法：**把约束转换成cost，这个cost由一个barrier function表示，在违反约束的时候给一个很大的cost。

如常用的$log$障碍函数

• 这是比较常用的方法，对于中小规模（变量范围在几千几万左右的）的凸优化问题来说这种方法是gold standard。
• 对于非凸问题，需要很多工程上的hacks和tricks来保证生效。

罚函数法

**核心想法：**把约束转成惩罚项，同样加到cost中。

• 实现很简单，但是难以达到高精度，无法保证最终的约束得到满足
• 会使整个KKT系统的Hessian矩阵变得病态ill-conditioned(由于这些非线性优化项引入的系统的Hessian，会使Hessian阵的特征值分布广，使条件数很大)，会让牛顿法的迭代（也不止牛顿法）变得更加慢。
• 要让约束完全得到满足，则约束的权重$\rho$要取到无穷大，这在实际情况是不能接受的，也会造成矩阵病态

ps：课上有同学提问，能否使用罚函数法或者内点法求解一段时间，再利用这个结果切换到积极点法中求解。老师回答说，这是工程上常用的技巧，先使用惩罚项来polish问题。

增广拉格朗日法（Augmented Lagrangian）

由于罚函数法存在的一些问题，我们转向增广拉格朗日法，核心想法是：我们在内层优化关于x的函数，在外层更新拉格朗日乘子来卸载（offloading）惩罚项

增广拉格朗日函数为：
$$\begin{array}{r}\underset{x}{\min }f(x)-{\tilde{\lambda }}^{T}c(x)+\rho /2[\min (0,c(x)){]}^2\end{array}\text{\hspace{1em}(11)}$$
将该函数对$x$求导，有
$$\begin{array}{r}\frac{\partial f}{\partial x}-\tilde{\lambda }\frac{\partial c}{\partial x}+\rho c(x{)}^T\frac{\partial c}{\partial x}=\frac{\partial f}{\partial x}-[\tilde{\lambda }-\rho c(x){]}^T\frac{\partial c}{\partial x}=0\\ ⟹\tilde{\lambda }\leftarrow \tilde{\lambda }-\rho c(x)\text{对于active的约束}.\end{array}$$

在初始时，令${\tilde{\lambda }}^T=0$，表示所有约束都不active，而后做一次无约束的最优化（关于$x$），判断有哪些约束被违反了（active），若有违反，则更新对应的拉格朗日乘子（拉格朗日乘子有值后才会在${\tilde{\lambda }}^{T}c(x)$项中存在惩罚），则在下一次的最优化中，就会有满足这个约束的趋势，一直迭代到收敛。

因为$c(x)\ge 0$是约束，在约束违反时候（$c(x)<0$），$\lambda$有一个正值，就会倾向于把$c(x)$优化地越大越好，这样就会有一定的裕度。

外层$\rho$的更新是为了不断增大惩罚力度，$\rho$ 越大则$\lambda$也会越来越大，最终当然取到无穷大的时候若存在解则不等式约束一定会被满足。当然我们期望的是不会取到无穷大，因为这样同样会使优化问题变得病态

案例分析

例子中是一个经典的QP问题，利用增广拉格朗日法解不等式约束。

function newton_solve(x0,λ,ρ)   # 这个函数不改变lambda和rho
    x = x0           # 问题是min x^2+2*y^2 约束是x-y+1=0 问题是二维的
    p = max.(0,c(x)) # 注意这个p和ρ是不同的，并且在这里例子中，约束是 c(x) <= 0
    C = zeros(1,2)   # 所以要cost中的λ是正号才能惩罚违反约束的情况（c(x) > 0）
    if c(x) ≥ 0      # 由于问题只有一个约束，若违反的话，那就设置梯度，用于优化
        C = ∂c(x)    
    end
    g = ∇f(x) + (λ+ρ*p)*C'  # 最终要优化的是关于x的拉格朗日函数，因此要对其求导，然后找根
    while norm(g) ≥ 1e-8
        H = ∇2f(x) + ρ*C'*C # Hessian阵，内层环不优化lambda，因此只有x
        Δx = -H\g
        x = x+Δx
        p = max.(0,c(x))    # 对应问题的最后一项
        C = zeros(1,2)
        if c(x) ≥ 0
            C = ∂c(x)       # 同样迭代完后判断有没有违反。没有违反就不要继续优化了，消掉COST和梯度
        end
        g = ∇f(x) + (λ+ρ*p)*C'
    end
    return x
end

增广拉格朗日法 初始点

增广拉格朗日法 迭代

在这个例子中，即使外层不更新$\rho$也能收敛

小结

• 相较于内点法/障碍函数法，这种方法不要求初始点一定要满足约束，即使初始点不满足约束，则后面会通过$\lambda$的增大逐渐把状态拉回来
• 可以很好地处理非凸问题，比较鲁棒，但是缺点就是没有特别快，但是也是超线性收敛的，
• 相较于罚函数法，这种方法可以用有限的$\rho$来解决问题，优化问题更好解（这点存疑，因为罚函数法好像也动态增加$\rho$）。另外一个优势是$\lambda$的更新是与$c(x)$的大小相关的，就是相当于有个自适应的过程，若$c(x)$的某些分量违反地多，那么对应地$\lambda$则大一些，其他违反地少地分量对应地$\lambda$则小一些。