数值优化（Numerical Optimization）(5)约束优化（三）

养生的控制人

发布时间 2022.03.28阅读数 3606 评论数 0

æ°å¼ä¼åï¼Numerical Optimizationï¼(5)çº¦æä¼åï¼ä¸ï¼

本次博客还是继续之前的约束优化问题，介绍求解线性规划的另一种方法——内点法（叫做内点法的原因是优化命题解的迭代过程是从约束域内部出发，沿着中心路径逐步走到边界），包含障碍函数法和原始对偶法。PS：实际应用时考虑求解效率和精度问题，通常首选的是原始对偶法。

参考自：

S. Boyd, L. Vandenberghe, Convex Optimization

2. CMU 10-725,Convex Optimization

3. 学弱猹的博客

障碍函数法

考虑问题

$\min f_0(x)\\ s.t. f_i(x) \le 0, i=1,...,m\\ Ax = b$

由于不等式约束的存在使得问题不能直接求解，那么就想办法去掉不等式约束，这里的做法是对优化问题进行改写使得优化命题不包含不等式约束，而目标函数中包含对不等式约束的惩罚，一个很直接的想法就是对不满足不等约束的情况进行惩罚，因此可以引入如下惩罚函数

${I_ - }(x) = 0,x \le 0\\ {I_ - }(x) =\infty ,x > 0$

当满足约束 $f_i(x)\le 0$ 时， $I_-(f_i(x))=0$ 即没有惩罚，当 $f_i(x)>0$ ，也就是说明我们不希望这种情况发生。利用惩罚函数代替不等式约束可以得到优化命题

$\min f_0(x) + \sum_{i=1}^m I_-(f_i(x))\\ s.t. Ax=b$

需要注意，这个目标函数是不连续的，为了使得目标函数连续，我们可以使用对数函数 $\log(-x)$ 来近似这个惩罚函数，因为对数函数在 $x$ 趋于 $0$ 时的函数值趋于无穷（惩罚无限大，也是我们想要的），同时还可以保证目标函数二次可微。定义 $\phi(x)=-\sum^{m}_{i=1} \log(-f_i(x))$ ，并引入系数 $1/t$ 来衡量近似程度，则上述优化命题转化为

$\min f_0(x) + \frac{1}{t}\phi(x)\\ s.t. Ax=b$

其中 $t$ 越大，用对数函数 $\log(-x)$ 来近似原来不可导的惩罚函数 $I_-(x)$ 的近似效果越好。

此时，原始带不等式约束的优化命题转化成为等式约束的优化命题，接下来我们来看看这一类问题怎么解。

将等数约束优化命题的目标函数乘以 $t$ ，（此时认为
$t$ 为固定的参数，不是变量）该问题的KKT条件可以表示为

$t\nabla f_0(x^*(t)) - \sum_{i=1}^m \frac{1}{f_i(x^*(t))} \nabla f_i(x^*(t)) + A^Tw = 0\\ Ax^*(t) = b\\ f_i(x^*(t)) \le 0, i=1,...,m$

可以证明给定 $x^*(t)$ 和对应的 $w$ ，（断句） $u_i^*(t)=-\frac{1}{tf_i(x^*(t))},i=1,...,m$ 和 $v^*(t)=w/t$ 是原问题的对偶可行解，其中原问题的拉格朗日函数为

$L(x,u,v) = f_0(x) + \sum_{i=1}^m u_if_i(x) + v^T(Ax-b)$

进一步可以证明对偶间隔满足

$f(x^*(t))-f^* \le m/t$

这可以作为算法的收敛准则，同时意味着随着 $t$ 趋于无穷， $x^*(t)$ 趋于 $x^*$ 。

因此求解上述KKT系统等价于求解方程组（其中 $v = w/t$ ）

$r(x,v) = \left[ {\begin{array}{*{20}{c}} {\nabla {f_0}(x) + \sum\limits_{i = 1}^m {( - \frac{1}{{t{f_i}(x)}})\nabla {f_i}(x)} + {A^T}v}\\ {Ax - b} \end{array}} \right] = 0$

根据一阶近似（ $y=(x,v)$ ）

$r(y+\Delta y) \approx r(y) + d r(y) \Delta y$

可以得到更新公式为（ $d$ 表示求导数）

$dr(y) \times \Delta y = -r(y)$

带入 $r(x,v)$ 的表达式，即左侧的导数对应为雅可比矩阵

$\left[ {\begin{array}{*{20}{c}} {{\nabla ^2}{f_0}(x) + \sum\limits_{i = 1}^m {\frac{1}{{r{f_i}{{(x)}^2}}}\nabla {f_i}(x)\nabla {f_i}{{(x)}^T}} - \sum\limits_{i = 1}^m {\frac{1}{{r{f_i}(x)}}{\nabla ^2}{f_i}(x)} }&{{A^T}}\\ A&0 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\Delta x}\\ {\Delta v} \end{array}} \right] = - r(x,v)$

这样就可以得到更新公式

$x_{k+1}=x_k + \Delta x \\ v_{k+1}=v_k + \Delta v$

每次更新解后还要增加 $t$ 的值，直到 duality gap 小于一定阈值即可。理论上来说一个 $t$ 可以对应一个解，假如说我们让 $t$ 连续变换，那么这个解也会连续变换，就会形成一条弧线。这一条弧线我们叫中心路径（central path）。

罚函数法的算法流程可以总结为（来自S. Boyd, L. Vandenberghe, Convex Optimization）

给定严格可行解 $x$ , $t:=t^0 >0$

重复一下步骤：

求解最小化问题 $\min tf_0+\phi$ s.t. $Ax=b$ 的最优解 $x^*(t)$
更新 $x:= x^*(t)$
如果 $m/t \le \epsilon$ 则停止，否则进入下一步
增加 $t:=\mu t$

原始对偶内点法

还是考虑刚才那个原始问题

$\min f_0(x)\\ s.t. f_i(x) \le 0, i=1,...,m \\ Ax = b$

该问题的扰动KKT条件（与原始KKT条件的不同之处在于互补松弛条件的等式右端是 $-\frac{1}{t}$ 而不是 $0$ ）为

$\nabla f_0(x)+\sum_{i=1}^m u_i \nabla f_i(x)+A^Tv =0\\ u_i\cdot f_i(x)=-\frac{1}{t},i=1,...,m\\ f_i(x) \le 0,i=1,...,m\\ Ax = b\\ u_i \ge 0$

定义

$F(x) = \left( {\begin{array}{*{20}{c}} {{f_1}(x)}\\ \vdots \\ {{f_m}(x)} \end{array}} \right),dF(x) = \left[ {\begin{array}{*{20}{c}} {\nabla {f_1}{{(x)}^T}}\\ \vdots \\ {\nabla {f_m}{{(x)}^T}} \end{array}} \right]$

KKT条件可以写成矩阵形式的方程组

$r(x,u,v)=\left( {\begin{array}{*{20}{c}} {\nabla {f_0}(x) + dF{{(x)}^T}u + {A^T}v}\\ { - diag(u)F(x) - \frac{1}{t} \cdot 1}\\ {Ax - b} \end{array}} \right) = 0$

这是一个非线性的系统 $r(x,u,v)=0$ ，仍然还是用牛顿求根的方法来计算，记 $y=\{x,u,v\}$ ，则 $r(y+\Delta y) \approx r(y) + d r(y) \Delta y$ （这里的 $d$ 代表求导），解的更新公式可以表示为

$dr(y) \times \Delta y = -r(y)$

定义符号

$r_{dual}=\nabla f_0(x) + dF(x)^Tu+A^Tv \\ r_{cent}=- diag(u)F(x) - \frac{1}{t} \cdot 1\\ r_{prim}=Ax-b$

表示点 $y=(x,u,v)$ 的对偶、中心和原始残差。解的更新公式表示为（左侧为雅可比矩阵）

$\left[ {\begin{array}{*{20}{c}} {{\nabla ^2}{f_0}(x) + \sum\limits_{i = 1}^m {{u_i}{\nabla ^2}{f_i}(x)} }&{DF{{(x)}^T}}&{{A^T}}\\ { - diag(u)DF(x)}&{ - diag(F(x))}&0\\ A&0&0 \end{array}} \right]\left( {\begin{array}{*{20}{c}} {\Delta x}\\ {\Delta u}\\ {\Delta v} \end{array}} \right) = - \left( {\begin{array}{*{20}{c}} {{r_{dual}}}\\ {{r_{cent}}}\\ {{r_{prim}}} \end{array}} \right)$

在罚函数法中，我们的对偶间隔为 $m/t$ ，而在原始对偶内点法中，我们构造一个替代对偶间隔

$\eta = -F(x)^Tu = -\sum_{i=1}^m u_i f_i(x)$

由于 $u_i = -1/(tf_i(x))$ 所以对应的 $\eta = m/t$

综上总结原始对偶内点法的步骤为（CMU/convexopt）

选择初始点 $x^{(0)}$ 使得 $f_i(x^{(0)}) < 0$ ， $u^{(0)} >0$ 。定义 $\eta^{(0)}=-F(x^{(0)})^Tu^{(0)}$ ，固定参数 $\mu >1$

重复以下步骤：

计算 $t=\mu m/\eta^{(k-1)}$
计算原始对偶更新方向 $\Delta y$
使用回溯法确定步长 $s$
更新解 $y^{(k)} = y^{(k-1)}+s \Delta y$
计算替代对偶间隔 $\eta^{(k)}=-F(x^{(k)})^T u^{(k)}$
判断替代对偶间隔是否小于一定阈值且原始、对偶残差的平方和 $(||r_{prim}||_2^2 + ||r_{dual}||_2^2)^{1/2}$ 小于一定阈值