【学界/编码】凸优化算法 I: 内点法(interior point method)求解线性规划问题

super1 · 发表于 2021-7-10 14:54

作者： @李恩志 Louisiana State University 计算物理学博士。

『运筹OR帷幄』责任编辑： @王源（东北大学系统工程博士生）

本篇文章是由以上作者在知乎上的优秀文章（原文链接:凸优化算法 I: 内点法(interior point method)求解线性规划问题），通过『运筹OR帷幄』责任编辑整理修改而成的。
欢迎原链接转发，付费转载请前往 @留德华叫兽 的主页获取信息，盗版必究。
敬请关注和扩散本专栏及同名公众号，会邀请全球知名学者陆续发布运筹学、人工智能中优化理论等相关干货、知乎Live及行业动态：
『运筹OR帷幄』大数据人工智能时代的运筹学

0前言

内点法是一种处理带约束优化问题的方法，其在线性规划，二次规划，非线性规划等问题上都有着很好的表现。在线性规划的问题上，相对于鼎鼎大名的单纯形法，内点法是多项式算法，而单纯形法并非多项式算法。从实际应用的效果来说，内点法也达到了足以和单纯形法分庭抗衡的地步，尤其针对大规模的线性规划问题内点法有着更大的发展潜力。

内点法并不仅仅用于线性规划的求解，值得一提的是内点法的很多思想有着更广泛的应用，例如障碍函数法的思想。

1导论

线性规划问题的一般形式为

$\text{min } c^{T}x \\ \text{s.t. } Ax \preceq b$
这里，目标函数为线性函数 $c^{T}x: \mathbb{R}^{n} \rightarrow \mathbb{R}$ ，约束条件为 $A_{ij}x_{j} \le b_{i}, i = 1, 2, ..., m.$ 矩阵  为 $m\times n$ 的满秩矩阵，其中 $m$ 为约束条件的个数， $n$ 为变量的个数。通常约束条件的个数大于变量的个数，所以有 $m > n$ .

2线性规划问题的等价(近似)表述

这个线性规划问题可以重新表述为计算 $\text{min } f(x)$ ，其中

$f(x) = c^{T}x + \sum_{i = 1}^{m}I(A_{ij}x_{j} - b_{i})$
这里，我们使用了一个indicator函数，定义为

$I(u) = \begin{cases} 0 & \text{if } x \le 0 \\ \infty & \text{if } x > 0 \end{cases}$
引入这个函数的意义在于可以将约束条件直接写入到目标函数里面，这样我们直接求新的函数的极小值就可以了，而不必借助于未知乘子。但是这里有一个问题，那就是indicator函数存在不可求导的点，因此在求函数极小值的时候我们没法通过普通的微分法来确定函数的极小值。为了规避这个问题，我们可以用一个光滑的函数来近似这个indicator函数。一个不错的选择是用 $I_{t}(u) = -\frac{1}{t}\log (-u)$ 来代替indicator函数。  只有在 $u < 0$ 的时候有定义，我们规定当 $u > 0$ 的时候 $I_{t}(u) = \infty$ . 而且参数 $t > 0$ 越大，函数  就越接近于 $I(u)$ . 所以我们可以通过调节  的值来调节这个函数的近似程度。使用这个近似的indicator函数，我们新的的目标函数可以写作

$f(x) =t c^{T}x - \sum_{i = 1}^{m} \log(-A_{ij} x_j + b_i)$ .
因为线性函数是凸函数，并且  也是凸函数，所以  是凸函数，因此我们可以很容易用凸优化的经典方法得到该函数的极小值。

3计算函数的梯度和Hessian矩阵

为了求函数的极小值，根据微积分的经典结果，我们只需令函数的梯度等于零，然后计算梯度为零时对应的解 $x^{\star}$ .
函数的梯度为

$\frac{\partial f}{\partial x_{k}} = t c_{k} + \sum_{i = 1}^{m}\frac{-A_{ik}}{A_{ij}x_{j} - b_{i}}$
Hessian 矩阵为

$\frac{\partial^2 f}{\partial x_{k}\partial x_{l}} = \sum_{i = 1}^{m}\frac{A_{ik} A_{il}}{\Big(A_{ij} x_{j} - b_{i} \Big)^2}$
定义对角型矩阵为

$D_{ij} = \delta_{ij} \frac{1}{\Big(A_{ik}x_{k} - b_{i}\Big)^2}$
于是Hessian矩阵可以写作

$H_{f} = A^{T} D A$
因为  为正定矩阵，所以Hessian矩阵至少为半正定矩阵。所以函数  是一个凸函数。而且矩阵  为可逆矩阵，矩阵  满秩，所以Hessian矩阵为可逆矩阵。于是函数  为强凸函数。所以，要计算  的根，我们可以用高效的牛顿迭代法。

4牛顿迭代法

我们现在的目标是计算  的根。因为Hessian矩阵可逆，所以我们可以用牛顿迭代法求解。牛顿迭代法为

$x^{(n+1)} = x^{(n)} - H^{-1}\nabla f$
在这个迭代过程中，参数  为固定的。每对应一个  ，我们都可以得到一个解  . 如果我们扫描参数  ，我们就可以得到一系列的解。其中最大的  的对应的解应该最精确。
一个更好的算法是选取一个比较小的初始参数  ，求出这个参数对应的解  . 然后增加  ，用之前得到的解  来初始化当前  所对应的牛顿迭代法的试解。这样算出来的解应该比直接计算  所对应的解更加精确。这样逐步迭代从小  到大  可以得到一系列的解，最后得到的解  称作是淬火解。这个解应该可以满足我们的精度需求。

5一个例子

为了展示该算法的功效，我们举一个简单的例子。这个例子可以用简单的几何方法求解出来。我们要展示的是，我们可以用数值方法得到同样的结果（在一定的精度范围之内）。之所以要用数值方法求解这个简单的问题，是因为数值方法对更复杂的问题同样有效，而简单的几何方法对复杂问题却已经不适用了。
现在要研究的例子为

$\text{min } x + y \\ \text{s.t. }\\ x + 2y \le 1 \\ 2x + y \le 1 \\ x \ge 0 \\ y \ge 0$
也就是

$c = \begin{pmatrix} 1 \\ 1 \end{pmatrix}$ , $A = \begin{pmatrix} 1 & 2 \\ 2 & 1 \\ -1 & 0 \\ 0 & -1 \end{pmatrix}$ , $b = \begin{pmatrix} 1 \\ 1 \\ 0 \\ 0 \end{pmatrix}$ .
梯度为

$\nabla f(x) = tc - A^{T} \begin{pmatrix} \frac{1}{x + 2y -1}\\ \frac{1}{2x + y -1} \\ \frac{1}{-x} \\ \frac{1}{-y} \end{pmatrix}$
Hessian矩阵为

$H_{f}(x) = A^{T}DA$
其中，对角型矩阵  为

$D = \begin{pmatrix} \frac{1}{(x + 2y - 1)^2 } & 0 & 0 & 0 \\ 0 & \frac{1}{(2x + y - 1)^2 } & 0 & 0 \\ 0& 0 & \frac{1}{x^2} & 0 \\ 0 & 0 & 0 & \frac{1}{y^2} \end{pmatrix}$
对于一个固定的参数  ，选择一个恰当的初始解 $x^{(0)}$ ，代入牛顿迭代公式

$x^{(n+1)} = x^{(n)} - H_{f}(x^{(n)})^{-1}\nabla f(x^{(n)}), n \ge 0$
可以得到一个依赖于参数  的解 $x^{\star}_{t}$ .
用几何方法很容易求得这个例子的解为 $(x^{\star}, y^{\star} ) = (0, 0)$ . 所以我们期待，当 $\lim_{t \rightarrow \infty} (x_{t}^{\star}, y_{t}^{\star}) = (0, 0)$ .
我需要用程序求这个例子的数值解。我写了一个Python程序实现这个算法，程序地址为PrimerLi/linear-programming.

程序结果如下：
这里用了四个不同的初始点来初始化程序，最终得到解都收敛到了 $(0, 0)$ . 箭头方向为参数  增加时 $(x_t^{\star}, y_t^{\star} )$ 移动的方向。这正是我们期待的结果。

在此感谢『运筹OR帷幄』审稿人对本文提出了宝贵的意见。
『运筹OR帷幄』审稿人 @覃含章，美国麻省理工学院（MIT）计算科学与工程方向博士在读，清华大学工业工程及数学与应用数学（第二学位）本科。研究兴趣主要为优化理论，机器学习算法在运营管理中的应用。
以上『运筹OR帷幄』专栏所有文章都会同步发送至留德华叫兽的头条主页，以及同名微信公众号，目前预计受众10w +
<hr/>如果你是运筹学/人工智能硕博或在读，请在下图的公众号后台留言：“加微信群”。系统会自动辨认你的关键字，并提示您进一步的加群要求和步骤，邀请您进全球运筹或AI学者群（群内学界、业界大佬云集）。
运筹学|控制论|优化理论爱好者，欢迎加qq群：686387574
人工智能爱好者，欢迎加qq群： 685839321
数据科学|分析爱好者，欢迎加qq群：130414574
最后敬请关注和扩散本专栏及同名公众号，会陆续发布运筹学、人工智能中优化理论相关干货及行业动态：『运筹OR帷幄』大数据和人工智能时代下的运筹学 - 知乎专栏

rustum · 发表于 2021-7-10 14:58

写的很好，已经打印学习，有四个问题请教一下，一是初始点是不是需要严格可行？第二个问题，每一次是不是把所有的不等式约束考虑，而不是和积极集等方法那样，只取有效不等式部分？第三，每一次在迭代过程中没有取到等号的不等式有没有什么数值上的特点？最后一个问题，迭代的过程中拉格朗日乘子能不能给出来啊？

JoshWindsor · 发表于 2021-7-10 14:59

加微信群

stonstad · 发表于 2021-7-10 15:00

谢谢

jquave · 发表于 2021-7-10 15:04

是在文末二维码的微信公众号后台回复哈

yukamu · 发表于 2021-7-10 15:04

加微信群

pc8888888 · 发表于 2021-7-10 15:04

这个线性规划问题可以重新表述为计算minf(x) 这里是否错误了，应该用max f(x) 才对吧

mypro334 · 发表于 2021-7-10 15:13

加微信群

unityloverz · 发表于 2021-7-10 15:13

牛顿迭代法中的H是指H_f吗

mypro334 · 发表于 2021-7-10 15:18

就是海森矩阵

		自动登录	找回密码
密码			立即注册

【学界/编码】凸优化算法 I: 内点法(interior point method)求解线性规划问题

本帖子中包含更多资源