线性规划-内点法初探

参考：最优化理论与算法（第2版）（陈宝林）

书中首先介绍了将一般线性规划问题转化为Karmarkar标准问题求解，为简化计算，Karmarkar等人又给出了内点法以求解线性规划问题，此部分在书中为*号引申内容，介绍较为简略，此处也是对书中内容做简要的补充。

考虑如下线性规划问题
$$\begin{array}{ll}\max & {\mathbf{c}}^{\mathrm{T}}\mathbf{x}\\ \text{ s. t. }& \mathbf{A}\mathbf{x}⩽\mathbf{b}\end{array}$$
其中：$\mathbf{c},x\in {\mathbb{R}}^n$，$\mathbf{A}\text{ 是 }m\times n\text{矩阵, }m⩾n\text{. }$
算法的基本思想,是从内点 ${\mathbf{x}}^{(0)}$ 出发,沿可行方向求出使目标函数值上升的后继点,再从得到的内点出发, 沿另一个可行方向求使目标函数值上升的内点。

将问题进行松弛：
$$\begin{array}{ll}\max & {\mathbf{c}}^{\mathrm{T}}\mathbf{x}\\ \text{ s. t. }& \mathbf{Ax}+\mathbf{v}=\mathbf{b},\\ & \mathbf{v}⩾0.\end{array}$$
对于第k轮迭代，即有
$${\mathbf{v}}^{(k)}=\mathbf{b}-\mathbf{A}{\mathbf{x}}^{(k)}$$

• 随后进行Affine Scaling，这是因为：当 ${\mathbf{x}}^k$ 是非常接近边界又不是最优解时, 步长将被迫选得非常小, 到最优解的收敛将非常慢。
• 故：若当前解 ${\mathbf{x}}^k$ 不是很靠近 “中心”, 需要将坐标重新拉伸 (re-scale), (仿射) 变换到靠近 “中心”的位置。
• 如何定义一个可行域的中心：
  若 ${\mathbf{x}}^k=\mathbf{e}$, 则
  (1) ${\mathbf{x}}^k$ 距离边界 1 个单位.
  (2) 因此只要步长 ${\alpha }^k<1$, 则可确保 ${\mathbf{x}}^{k+1}>0$.

故可以定义对角矩阵：${\mathbf{D}}_k=\mathrm{diag}\left(\frac{1}{v_1^{(k)}},\cdots ,\frac{1}{v_m^{(k)}}\right).$，利用其进行Affine Scaling：$\mathbf{w}={\mathbf{D}}_k\mathbf{v}$，原问题可进一步被写为：
$$\begin{array}{ll}\max & {\mathbf{c}}^{\top }\mathbf{x}\\ \text{ s. t. }& \mathbf{Ax}+{\mathbf{D}}_k^{-1}\mathbf{w}=\mathbf{b},\\ & \mathbf{w}⩾0.\end{array}$$

在变换后的空间中，定义搜索方向为：
$$d=\left[\begin{array}{l}{d}_x\\ d_w\end{array}\right]\in {\mathbb{R}}^{m+n}$$
代入可得：
$$\begin{array}{rl}& A(x+d_x)+D_k^{-1}\left(w+d_w\right)=b\\ & \underset{b}{\underbrace{Ax+D_k^{-1}w}}+A{d}_x+D_k^{-1}{d}_w=b\\ & D_{k}A{d}_x+d_w=0\end{array}$$
随后，为了寻找一个“好方向”，
对于满足$D_{k}A{d}_x+d_w=0$的任一解，有 $A^{\mathrm{T}}{D}_k\left(D_{k}A{d}_x+d_w\right)=0$ (零空间投影)
可以得到：
$$d_x=-{\left(A^{\top }{D}_k^{2}A\right)}^{-1}{A}^{\top }{D}_k{d}_w.$$
为了使$c^T(x+d_x)$最大，将上式代入$c^T{d}_x$得到：
$${\mathbf{c}}^{\mathrm{T}}{\mathbf{d}}_x={\mathbf{c}}^{\mathrm{T}}\left[-{\left({\mathbf{A}}^{\mathrm{T}}{\mathbf{D}}_k^2\mathbf{A}\right)}^{-1}{\mathbf{A}}^{\mathrm{T}}{\mathbf{D}}_k{\mathbf{d}}_w\right]=-{\left[{\mathbf{D}}_k\mathbf{A}{\left({\mathbf{A}}^{\mathrm{T}}{\mathbf{D}}_k^2\mathbf{A}\right)}^{-1}\mathbf{c}\right]}^{\mathrm{T}}{\mathbf{d}}_w.$$
为了使上式最大，需进一步确定$d_w$取值，即选择同方向：
$${\mathbf{d}}_{\mathrm{w}}=-{\mathbf{D}}_k\mathbf{A}{\left({\mathbf{A}}^{\mathrm{T}}{\mathbf{D}}_k^2\mathbf{A}\right)}^{-1}\mathbf{c}.$$
又由于：$D_{k}A{d}_x+d_w=0$ 可知
$$-D_{k}A{d}_x=d_w=-D_k\mathbf{A}{\left(A^{\mathrm{T}}{D}_k^{2}A\right)}^{-1}c$$
可得到$$d_x={\left({\mathbf{A}}^{\mathrm{T}}{\mathbf{D}}_k^2\mathbf{A}\right)}^{-1}\mathbf{c}$$
随后，将得到的$d_w$逆仿射：
$$d_v=D_k^{-1}{d}_w=-A{\left(A^{\mathrm{T}}{D}_k^{2}A\right)}^{-1}c=-A{d}_x$$
故迭代更新表示为：
$${\mathbf{x}}^{(k+1)}={\mathbf{x}}^{(k)}+\lambda {\mathbf{d}}_{\mathbf{x}}$$
为保证更新序列始终为内点：
$$\begin{array}{rl}& \mathbf{A}\left({\mathbf{x}}^{(k)}+\lambda {\mathbf{d}}_x\right)<\mathbf{b},\\ & \lambda \mathbf{A}{\mathbf{d}}_x<\mathbf{b}-\mathbf{A}{\mathbf{x}}^{(k)}\\ & -\lambda {\mathbf{d}}_v<{\mathbf{v}}^{(k)}\end{array}$$
令：
α = min ⁡ { v i ( k ) − ( d v ) i ∣   ( d v ) i < 0 , i ∈ { 1 , ⋯   , m } } , \alpha=\min \left\{\left.\frac{v_i^{(k)}}{-\left(d_v\right)_i} \right\rvert\,\left(d_v\right)_i<0, i \in\{1, \cdots, m\}\right\}, α=min{−(dv​)i​vi(k)​​ ​(dv​)i​<0,i∈{1,⋯,m}},
取 $$\lambda =\gamma \alpha ，\gamma \in (0,1)$$
这样即可从$x^{(k)}$出发沿方向$d_x$求得使$c^{T}x$上升的点