内点法初探——线性规划标准形式下的求解思路

一般的线性规划具有以下形式：

\begin{alignat}{2} \min_{x} \quad & c^Tx + d\\ \mbox{s.t.}\quad &Gx \leq h\\ &Ax=b\\ \end{alignat}

其中，线性规划标准形是线性规划的一种特殊情况，近年来已经被广泛、深入地研究。在求解线性规划问题时，可以将上述的一般形式通过某种变化（如引入松弛变量等）转换成标准形式：

\begin{alignat}{2} \min_{x} \quad & c^Tx\\ \mbox{s.t.}\quad &Ax=b\\ &x \geq 0\\ \end{alignat}

其中

x\in R^n,A\in R^{m\times n},b\in R^m

本文主要讨论利用内点法求解线性规划标准形的过程。

内点法求解线性等式和不等式约束的优化问题，是通过将其简化成一系列线性等式约束问题求解。

首先，重新表述标准形问题，把不等式约束隐含在目标函数中：

\begin{alignat}{2} \min_{x} \quad & c^Tx + \sum_{i=1}^nI_{-}(-x_i)\\ \mbox{s.t.}\quad &Ax=b\\ \end{alignat}

其中Indicator函数

I_{-}(u)=\begin{equation} \left\{ \begin{array}{**lr**} 0 & u\leq 0\\ \infty, & u > 0 \end{array} \right. \end{equation}

不可微，因此需要查找一个替代函数来近似Indicator函数。

障碍法（Barrier Method）

障碍法的基本思想是用以下的函数近似Indicator函数：

I_t(u)=-\frac{1}{t}log(-u)

t越大越逼近Indicator函数。

代入可得（为了方便，我们用t乘目标函数考虑等价问题）

\begin{alignat}{2} \min_{x} \quad & tc^Tx + \sum_{i=1}^n-log(x_i)\\ \mbox{s.t.}\quad &Ax=b \end{alignat}

在上述目标中引入Lagrange乘子构建对偶问题有：

L(x,v)= tc^Tx - \sum_{i=1}^nlog(x_i) + A^Tv

对应的KKT条件为

\begin{alignat}{2} tc-diag(x)^{-1}1+A^Tv=0 \quad (1)\\ Ax=b \quad (2) \end{alignat}

利用Newton Step可以有

[\begin{matrix} diag(x)^{-2} & A^T \\ A & 0 \end{matrix} ] [\begin{matrix} \Delta x \\ \Delta v\end{matrix}] = [\begin{matrix} -tc+diag(x)^{-1}1-A^Tv \\ 0\end{matrix}]

整理可得

[\begin{matrix} diag(x)^{-2} & A^T \\ A & 0 \end{matrix} ] [\begin{matrix} \Delta x \\ w\end{matrix}] = [\begin{matrix} -tc+diag(x)^{-1}1 \\ 0\end{matrix}]

其中

w=v+\Delta v

.

通常通过消去

\Delta x

来求解方程，从第一个等式可得

\Delta x=-tdiag(x)^2c+x-diag(x)^2A^Tw

带入第二个方程得

Adiag(x)^2A^Tw=-tAdiag(x)^2c+b

综上，使用barrier method求解标准形的线性规划问题的步骤可以整理如下：

step1: 初始化

t_0

和可行点

x_0

step2: 根据

Adiag(x)^2A^Tw=-tAdiag(x)^2c+b

求解

w

step3: 根据

\Delta x=-tdiag(x)^2c+x-diag(x)^2A^Tw

求解

\Delta x

step4: 更新

x=x+\Delta x

step5: 判断

m/t<\epsilon

是否成立，若成立，则退出循环，输出

x

，否则执行step6

step6: 更新

t=\mu t

，执行step2

原对偶内点法（Primal-dual interior-point method）

原对偶内点法对应的KKT条件与Barrier method方法类似：

\begin{alignat}{2} c+diag(-x)^{-1}u+A^Tv=0 \quad (1. r_{prim})\\ diag(u)(-x)=-\frac{1}{t}1\quad(2.r_{cent}) \\ Ax=b \quad (3.r_{dual}) \end{alignat}

稍微整理一下可得

\begin{alignat}{2} c-diag(x)^{-1}u+A^Tv=0 \quad (1)\\ diag(u)(x)=\frac{1}{t}1\quad(2) \\ Ax=b \quad (3) \end{alignat}

利用Newton Step可得

[\begin{matrix} 0 & -1 &A^T \\ diag(u) & diag(x) & 0 \\ A & 0 & 0 \end{matrix} ] [\begin{matrix} \Delta x \\ \Delta u \\ \Delta v\end{matrix}] = [\begin{matrix} -c+diag(x)^{-1}u-A^Tv \\ -diag(u)x+\frac{1}{t}1 \\ b-Ax\end{matrix}]

代入

\Delta u=-diag(x)^{-1}(diag(u)\Delta x+diag(u)x-\frac{1}{t}1)=-diag(\frac{u}{x})\Delta x-u+\frac{1}{t}\frac{1}{x}

可得

[\begin{matrix} diag(\frac{u}{x}) &A^T \\ A & 0 \end{matrix} ] [\begin{matrix} \Delta x \\ \Delta v\end{matrix}] = [\begin{matrix} -c+\frac{1}{t} diag(x)^{-1}1-A^Tv \\ b-Ax\end{matrix}]

进一步代入可得

\Delta x=-diag(\frac{x}{u})c+\frac{1}{t}diag(\frac{1}{u})1-diag(\frac{x}{u})A^Tv-diag(\frac{x}{u})A^T\Delta v

Adiag(\frac{x}{u})A^T\Delta v=Ax-b-Adiag(\frac{x}{u})c+\frac{1}{t}A\frac{1}{u}-Adiag(\frac{x}{u})A^Tv

依次计算

\Delta v,\Delta x,\Delta u

的值。

综上，使用primal dual求解标准形的线性规划问题的步骤可以整理如下：

step1: 初始化

u_0>0,v_0,x_0

，定义

\eta_0=x_0^Tu_0

，定义参数

\mu > 1

step2: 定义

t=\frac{\mu m}{\eta_{k-1}}

，计算

\Delta x, \Delta u,\Delta v

step3: 确定步长s，更新

x_k,u_k,v_k

step4: 计算

\eta_k=x_k^Tu_k

，判断

\eta_k \leq \delta

且

(||r_{prim}||_2^2+||r_{dual}||_2^2)^{\frac{1}{2}} \leq \delta

,退出循环，同时输出

x

，否则重复step2

齐次内点法（Homogeneous interior-point method）

这里主要介绍Mosek方法。

原问题的对偶问题可以表示为

\begin{alignat}{2} \max_{y} \quad & b^Ty\\ \mbox{s.t.}\quad &A^Ty+z=c\\ &s \geq 0\\ \end{alignat}

原问题的最优性条件表示为

\begin{alignat}{2} Ax=b,x \geq 0\quad (1)\\ A^Ty+z=c, z \geq 0 \quad (2) \\ diag(x)z=0 \quad (3) \end{alignat}

原问题和对偶问题的对偶间隔为

c^Tx-b^Ty=x^Tz \geq 0

引入两个非负变量

\tau,\kappa

，化简齐次模型得到HLF模型

\begin{alignat}{2} Ax-b\tau=0 \quad (1)\\ A^Ty+z-c\tau = 0\quad (2) \\ -c^Tx+b^Ty-\kappa=0 \quad (3) \\ x, \tau, z,\kappa \geq 0 \quad \quad \quad \end{alignat}

显然，0解是一个合理但是没什么用的解。

如果我们通过HLF模型得到一组解

(x^*,\tau^*,y^*,z^*,\kappa^*)

，那么原问题和对偶问题的解为

(x,y,z)=(x^*/\tau^*,y^*/\tau^*,z^*/\tau^*)

对偶间隔为

\kappa^*/\tau^*

。

求解HLF模型需要满足以下5个条件：

\begin{alignat}{2} Ax-b\tau=0 \quad (1)\\ A^Ty+z-c\tau = 0\quad (2) \\ -c^Tx+b^Ty-\kappa=0 \quad (3) \\ diag(x)z=\mu 1 \quad (4)\\ \tau \kappa = \mu \quad (5) \\ x, \tau, z,\kappa \geq 0 \quad \quad \quad \end{alignat}

对应残差为

\begin{alignat}{2} r_p=b\tau-Ax \quad (1)\\ r_d=c\tau - A^Ty-z \quad (2) \\ r_g=c^Tx-b^Ty+\kappa \quad (3) \\ \mu=(x;\tau)^T(z;\kappa)/(n+1) \quad (4) \end{alignat}

搜索更新方向为

\begin{alignat}{2} Ad_x-bd_\tau=\hat r_p=(1-\gamma)r_p \quad (1)\\ A^Td_y+d_z-cd_\tau = \hat r_d =(1-\gamma)r_d\quad (2) \\ -c^Td_x+b^Td_y-d_\kappa=\hat r_g=(1-\gamma)r_g \quad (3) \\ diag(s)d_x+diag(x)d_z=\hat r_{x,s}=-diag(x)z+\gamma \mu 1 \quad (4)\\ \kappa d_\tau + \tau d_\kappa = \hat r_{\tau,\kappa}=-\tau \kappa + \gamma\mu \quad (5) \end{alignat}

写成方程组形式

[\begin{matrix} A & -b & 0 & 0 & 0 \\ 0 & -c & A^T & I & 0 \\ -c^T & 0 & b^T & 0 & -1 \\ diag(z) & 0 & 0 & diag(x) & 0 \\ 0 & \kappa & 0 & 0 & \tau \end{matrix} ] [ \begin{matrix} d_x \\ d_\tau \\ d_y \\ d_z \\ d\kappa\end{matrix} ] =[\begin{matrix} \hat r_p \\ \hat r_d \\ \hat r_g \\ \hat r_{x,z} \\ \hat r_{\tau, \kappa }\end{matrix}]

代入

d_z=diag(x)^{-1}(\hat r_{x,z} - diag(z)d_x)

和

d_\kappa = \tau^{-1}(\hat r_{\tau,\kappa}-\kappa d_{\tau})

得

[\begin{matrix} -diag(s/x) & A^T & -c \\ A & 0 & b \\ -c^T & b^T & \kappa / \tau \end{matrix} ] [ \begin{matrix} d_x \\ d_y \\ d_\tau \end{matrix} ] =[\begin{matrix} \hat r_d-diag(x)^{-1}\hat r_{x,s} \\ \hat r_p \\ \hat r_g +\hat r_{\tau, \kappa } / \tau\end{matrix}]

定义

K=[\begin{matrix} -diag(z/x) & A^T \\ A & 0\end{matrix}]

通过求解

K(p;q)=(c;b)

和

K(u;v)=[\begin{matrix} \hat r_d - diag(x)^{-1} \hat r_{x,z} \\ \hat r_p \end{matrix}]

来计算

\begin{alignat}{2} d_\tau = \frac{\hat r_g + \hat r_{\tau,\kappa}/\tau - (-c;b)^T(u;v)}{\kappa/\tau + (-c;b)^T(p;q)} \quad (1) \\ d_x = u + p d_\tau \quad (2) \\ d_y = v + q d_\tau \quad (3) \end{alignat}

综上，使用mosek求解标准形的线性规划问题的步骤可以整理如下：

step1: 初始化

(x^0,\tau^0 , y^0, z^0,\kappa^0)=(e,1,0,e,1)

for k = 1,2, ...

step2: 计算

r_p^{k-1},r_d^{k-1},r_g^{k-1},\mu^{k-1}

，如果

max(\frac{||r_p^{k-1}||}{max(1,||r_p^0||)},\frac{||r_d^{k-1}||}{max(1,||r_d^0||)},\frac{|r_g^{k-1}|}{max(1,|r_g^0|)},\frac{|c^Tx^k/\tau^k-b^Ty^k/\tau^k|}{1+|b^Ty^k/\tau^k|} )<threshold

则输出

x^{k-1}/\tau^{k-1}

，如果

max(\frac{||r_p^{k-1}||}{max(1,||r_p^0||)},\frac{||r_d^{k-1}||}{max(1,||r_d^0||)},\frac{|r_g^{k-1}|}{max(1,|r_g^0|)})<threshold

且

\tau^{k-1}<threshold \times max(1,\kappa^{k-1})

，则表示该线性规划无解，退出循环。

step3: 初始化

\gamma = 0

，计算

\hat r_p,\hat r_d,\hat r_g, \hat r_{x,z}, \hat r_{\tau,\kappa}

step4: 更新

\gamma = (1-\alpha)^2min(\beta,1-\alpha)

，再次计算

\hat r_p,\hat r_d,\hat r_g, \hat r_{x,z}, \hat r_{\tau,\kappa}

step5: 更新

x(k),\tau(k) , y(k), z(k),\kappa(k)

，回step2

注：其中step4中的

\alpha

为搜索方向。