线性规划：内点法

本文介绍求解线性规划问题的内点法。它是一个多项式时间算法，在实际应用中效率也很高。尤其是对求解大规模线性规划，一些经验说，内点法比单纯形法更快。此外，内点法还可以被扩展，用来求解凸优化以及非线性规划问题。

考虑线性规划标准问题及其对偶问题：

原始问题（P）
$$\begin{array}{rl}\min & c^{T}x\\ \text{s.t. }& Ax=b\\ & x\ge 0\end{array}$$
对偶问题（D）
$$\begin{array}{rl}\max & b^{T}y\\ \text{s.t. }& A^{T}y+s=c\\ & s\ge 0\end{array}$$
其中 $A\in {\mathbb{R}}^{m\times n}$，$c,x,s\in {\mathbb{R}}^n$，$b,y\in {\mathbb{R}}^m$，且矩阵 $A$ 行满秩。

内点

先定义原始问题和对偶问题的可行域：
F ( P ) = { x ∈ R n : A x = b , x ≥ 0 } F ( D ) = { ( y , s ) ∈ R m × R n : A T y + s = c , s ≥ 0 } \begin{aligned} & \mathcal{F}(P) =\{x\in\mathbb{R}^n: Ax=b,x\geq0\}\\ & \mathcal{F}(D) = \{(y,s)\in\mathbb{R}^m\times\mathbb{R}^n: A^Ty+s=c, s\geq 0\} \end{aligned} ?F(P)={ x∈Rn:Ax=b,x≥0}F(D)={ (y,s)∈Rm×Rn:ATy+s=c,s≥0}?
接下来定义可行域的 内部（interior）：
F ? ( P ) = { x ∈ R n : A x = b , x > 0 } F ? ( D ) = { ( y , s ) ∈ R m × R n : A T y + s = c , s > 0 } \begin{aligned} & \mathcal{F}^{\circ}(P) =\{x\in\mathbb{R}^n: Ax=b, x > 0\}\\ & \mathcal{F}^{\circ}(D) = \{(y,s)\in\mathbb{R}^m\times\mathbb{R}^n: A^Ty+s=c, s > 0\} \end{aligned} ?F?(P)={ x∈Rn:Ax=b,x>0}F?(D)={ (y,s)∈Rm×Rn:ATy+s=c,s>0}?
本文介绍原始对偶内点法（Primal Dual Interior Point Method），它的思路是在可行域的内部 ${\mathcal{F}}^{?}(P)\times {\mathcal{F}}^{?}(D)$ 进行迭代，朝着目标优化的方向去逼近最优解。

在可行域内部迭代有什么好处？

给定内部的任意一点，存在它的一个 邻域。如果目标函数是光滑的，那么在内点可导，于是可以用牛顿法解方程。而这个待解的方程就是最优解的充分必要条件（见下文）。

罚函数

给定一个内点，我们想要迭代始终发生在可行域的内部。思路是引入 罚函数（也称障碍函数）：
$$F(x)=?\sum _{j=1}^n\ln (x_j)$$
给定 $\mu >0$，定义目标函数：
$$B_{\mu }(x)=c^{T}x+\mu F(x)$$
当 $x\to 0$ 时，$F(x)\to \infty$，换句话说，$B_{\mu }(x)$ 的定义域保证了 $x>0$。

接下来考虑在 ${\mathcal{F}}^{?}(P)$ 中最小化 $B_{\mu }(x)$，即
$$\begin{array}{rl}\min & c^{T}x+\mu F(x)\\ \text{s.t. }& Ax=b\\ & x>0\end{array}(\text{PP})$$
当 $\mu \to 0$ 时，$B(\mu )\to c^{T}x$，相当于最小化 $c^{T}x$；当 $\mu \to \infty$ 时，相当于最小化 $F(x)$。

下面介绍最优解的充要条件。

最优条件

设 ${\mathcal{F}}^{?}(P)$ 和 ${\mathcal{F}}^{?}(D)$ 非空。如果 $x$ 是上述问题 $(\text{PP})$ 的最优解，那么当且仅当存在 $(y,s)\in {\mathcal{F}}^{?}(D)$，满足如下条件：

1. 原始可行：$Ax=b$

2. 对偶可行：$A^{T}y+s=c$

3. 互补松弛：$XSe=\mu e$

其中 $X=diag(x)$，$S=diag(s)$，$e=(1,1,\dots ,1{)}^T\in {\mathbb{R}}^n$。

下面用等式描述上面的条件。

定义函数
$$F(x,y,s)=?\begin{array}{c}Ax?b\\ A^{T}y+s?c\\ XSe?\mu e\end{array}?$$
$x$ 是 $\text{(PP)}$ 的最优解等价于下面的方程有解：
$$F(x,y,s)=0$$

牛顿法

可以用牛顿法迭代法求解上面的方程 $F(x,y,s)=0$。

这里只介绍牛顿法的基本步骤。记 $x^k,y^k,s^k$ 代表第 $k$ 步的迭代结果，那么
$$(x^{k+1},y^{k+1},s^{k+1})=(x^k,y^k,s^k)+{\alpha }^k?(\Delta x^k,\Delta y^k,\Delta s^k)$$
其中 $(\Delta x^k,\Delta y^k,\Delta s^k)$ 代表迭代方向，${\alpha }^k\in (0,1]$ 代表步长。

令 $J(x,y,s)$ 代表 $F(x,y,s)$ 的雅可比矩阵，即
$$J(x,y,s)=?\begin{array}{ccc}\frac{?F_1}{?x}& \frac{?F_1}{?y}& \frac{?F_1}{?s}\\ \frac{?F_2}{?x}& \frac{?F_2}{?y}& \frac{?F_2}{?s}\\ \frac{?F_3}{?x}& \frac{?F_3}{?y}& \frac{?F_3}{?s}\end{array}?=?\begin{array}{ccc}A& 0& 0\\ 0& A^T& I\\ A& 0& 0\\ S& 0& X\end{array}?$$
求解下面的线性方程组：
$$J(x,y,s)?\begin{array}{c}\Delta x\\ \Delta y\\ \Delta s\end{array}?+F(x,y,s)=0$$
即，
$$?\begin{array}{ccc}A& 0& 0\\ 0& A^T& I\\ S& 0& X\end{array}??\begin{array}{c}\Delta x\\ \Delta y\\ \Delta s\end{array}?=??\begin{array}{c}Ax?b\\ A^{T}y+s?c\\ XSe?\mu e\end{array}?$$
可以得到牛顿法的迭代方向。

由于 $x,y$ 是可行解，上式可以写成
$$?\begin{array}{ccc}A& 0& 0\\ 0& A^T& I\\ S& 0& X\end{array}??\begin{array}{c}\Delta x\\ \Delta y\\ \Delta s\end{array}?=??\begin{array}{c}0\\ 0\\ XSe?\mu e\end{array}?$$

中心路径

回顾问题 $\text{(PP)}$，它在 ${\mathcal{F}}^{?}(P)$ 中最小化最小化 $B_{\mu }(x)$，即
$$\begin{array}{rl}\min & c^{T}x+\mu F(x)\\ \text{s.t. }& Ax=b\\ & x>0\end{array}(\text{PP})$$
我们已知它最优解的充要条件，而且可以用牛顿法求它的最优解。

但是，给定 $\mu >0$，$\text{(PP)}$ 的最优解并不是原问题 $\text{(P)}$ 的最优解。因此，迭代的目标不是求解 $\text{(PP)}$ ，而是去逼近原问题 $\text{(P)}$ 的最优解。逼近的方式就是在迭代中不断减小 $\mu$ 值，然后更新下降方向；当 $\mu$ 足够小时，得到的解与原问题的最优解足够接近。

具体来说，给定初始值 ${\mu }^0>0$ 和初始点 $x^0({\mu }^0),y^0({\mu }^0),s^0({\mu }^0)\in {\mathcal{F}}^{?}(P)\times {\mathcal{F}}^{?}(D)$ ，迭代过程得到如下的点列：
$$\begin{gathered} x^0({\mu }^0),y^0({\mu }^0),s^0({\mu }^0), \\ {x}^1({\mu }^1),y^1({\mu }^1),s^1({\mu }^1), \\ ? \\ {x}^k({\mu }^k),y^k({\mu }^k),s^k({\mu }^k) \end{gathered}$$
其中 ${\mu }^0>{\mu }^1>?>{\mu }^k$。当 ${\mu }^k<?$ 时，迭代停止，其中 $?>0$ 代表逼近的精度。

我们把 { x ( μ ) , y ( μ ) , s ( μ ) : μ > 0 } \{x(\mu),y(\mu),s(\mu): \mu > 0\} { x(μ),y(μ),s(μ):μ>0} 称为 中心路径（Central Path）。迭代的过程就是沿着中心路径，令 $\mu$ 值不断减小的过程，这样的算法也称为 路径跟踪（Path Following）。

内点法（理论版）

接下来介绍具体的迭代步骤。

原始对偶内点法 (Primal-Dual Interior Point Method)

---

第一步，初始化 ${\mu }^0,x^0,y^0,s^0$。注意： $\mu ,x,s>0$ 且 $x^0,y^0,s^0$ 可行。

第二步，判断误差：如果 $\mu <?$ 则停止。

第三步，计算牛顿方向：
$$?\begin{array}{ccc}A& 0& 0\\ 0& A^T& I\\ S& 0& X\end{array}??\begin{array}{c}\Delta x\\ \Delta y\\ \Delta s\end{array}?=??\begin{array}{c}0\\ 0\\ XSe?\mu e\end{array}?$$
第四步，更新 $x,y,s$：
$$(x,y,s)\leftarrow (x,y,s)+(\Delta x,\Delta y,\Delta s)$$
第五步，减小 $\mu$ 值，令 $\mu \leftarrow \sigma ?\mu$，其中 $\sigma \in [0,1]$， 然后执行第二步。

---

这里还有几个问题需要回答：

1. $\mu$ 和 $\sigma$ 值如何选取？

2. 第四步更新 $x,y,s$ 之后，它还是可行解吗？

3. 如何计算初始可行解 $x^0,y^0,s^0$ ？

先看第一个问题。

先计算原问题 $\text{(P)}$ 和对偶问题 $\text{(D)}$ 的对偶间隙 (Duality Gap)：
$$c^{T}x?b^{T}y=(A^{T}y+s{)}^{T}x?b^{T}y=s^{T}x$$
令 $\mu =\frac{s^{T}x}{n}$，当 $\mu <?$ 时，我们有 $c^{T}x?b^{T}y\le n?$，即目标函数值与最优值当误差不超过 $n?$。

此外，令 $\sigma =1?\frac{0.4}{\sqrt{n}}$，可以证明最多 $O(\sqrt{n}\ln (1/?))$ 次迭代，就可以使得 $\mu <?$。

再看第二个问题。

只需要验证 $x+\Delta x$ 与 $(y+\Delta y,s+\Delta s)$ 是否满足原始问题和对偶问题的约束：
$$\begin{array}{rl}& A(x+\Delta x)=b+A\Delta x=b\\ & A^T(y+\Delta y)+(s+\Delta s)=c+A^T\Delta y+\Delta s=c\end{array}$$
因此答案是肯定的。

关于第三个问题，计算初始可行解并不容易，这里不做介绍。原因是在实际应用中，可以通过逼近的方式满足可行性，因而初始解只要保证 $x,s>0$ 即可。

内点法（实践版）

编程实现内点法时，一般不会严格按照上面的理论版，原因是效率不高。

下面介绍一个实践版本。

---

第一步，初始化 $x^0>0,s^0>0,y^0$，令 $\mu =\frac{s^{T}x}{n}$。

第二步，判断误差：如果 $\mu <?$ 则停止，其中
μ = max ? { ∣ ∣ A x ? b ∣ ∣ , ∣ ∣ A T y + s ? c ∣ ∣ , ∣ ∣ s T x ∣ ∣ } \mu = \max \{||Ax-b||, ||A^Ty+s-c||, ||s^Tx||\} μ=max{ ∣∣Ax?b∣∣,∣∣ATy+s?c∣∣,∣∣sTx∣∣}
第三步，计算牛顿方向：
$$?\begin{array}{ccc}A& 0& 0\\ 0& A^T& I\\ S& 0& X\end{array}??\begin{array}{c}\Delta x\\ \Delta y\\ \Delta s\end{array}?=??\begin{array}{c}Ax?b\\ A^{T}y+s?c\\ XSe?\mu e\end{array}?$$
第四步，更新 $x,y,s$：
$$\begin{array}{rl}& x\leftarrow x+{\alpha }_P?\Delta x\\ & (y,s)\leftarrow (y,s)+{\alpha }_D?(\Delta y,\Delta s)\end{array}$$
其中 ${\alpha }_P,{\alpha }_D$ 使得 $x,s>0$，即
$$\begin{gathered} {\alpha }_P=\min \left\{1,\underset{\Delta x_j<0}{\min }\left\{\frac{x_j}{?\Delta x_j}\right\}\right\} \\ {\alpha }_D=\min \left\{1,\underset{\Delta s_j<0}{\min }\left\{\frac{s_j}{?\Delta s_j}\right\}\right\} \end{gathered}$$
第五步，减小 $\mu$ 值，令 $\mu \leftarrow \sigma ?\mu$，其中 $\sigma =0.1$ (经验值)， 然后执行第二步。

---

注意两点变化：

第一，初始解不要求可行，原因是可以通过迭代的方式保障可行性。它通过 $\mu$ 的选择来实现，它是三种误差的最大值，即原问题的可行性误差 $||Ax?b||$， 对偶问题的可行性误差 $||A^{T}y+s?c||$，以及目标函数误差 $s^{T}x$。 当 $\mu <?$ 且 $?$ 精度足够小时，原问题和对偶问题可行，并且目标函数值达到最优。

第二，由于迭代过程中，原始解和对偶解可能不满足可行性，因此更新 $x,y,s$ 的步长 $\alpha$ 需要保证 $x,s>0$。其中 $x$ 与 $(y,s)$ 也可以采用相同的步长，例如 $\alpha =\min ({\alpha }_D,{\alpha }_P)$。

算法实现

下面用 Python 实现内点法。

先定义问题的输入和输出。

import numpy as npclass IPMlp(object):""" A Primal Dual Interior Point Method for linear programming.说明：1. 最小化问题2. 行满秩"""def __init__(self, c, A, b):""":param c: n * 1 vector:param A: m * n matrix:param b: m * 1 vector"""# 输入self._c = np.array(c)self._A = np.array(A)self._b = np.array(b)# 输出self._obj = None  # objective function valueself._x = None  # primal solutionself._y = None  # dual solution# 辅助变量self._m = len(self._b)self._n = len(self._c)self._s = None  # dual slack variableself._alpha = Noneself._sigma = Noneself._mu = Noneself._iter_num = 0self._status = None

接下来是算法实现的思路。

class IPMlp(object):# 其它函数省略# ...def solve(self):self._init()  # initialization# 判断停止条件while not self._is_stop_criterion_satisfied():# 解线性方程组，得到牛顿下降方向dx, dy, ds = self._solve_linear_system()if dx is None:self._status = 'UNBOUNDED or INFEASIBLE'break# 计算步长，保证 x, s > 0self._alpha = self._get_alpha(dx, ds)# 更新原始解和对偶解# 这里采用了统一的步长self._x += self._alpha * dxself._y += self._alpha * dyself._s += self._alpha * ds# 减小 mu 值self._mu = (self._x @ self._s / self._n) * self._sigma# 更新目标函数self._obj = self._c @ self._x

更多细节请参考源码：完整代码

参考文献

[1]. David P. Williamson. ORIE 6300 Mathematical Programming I (lecture notes). 2014.

[2]. R. M. Freund and J. Vera. Interior-Point Methods for Linear Optimization. 2014.