Numerical Optimization Ch14. Linear Programming: Interior-Point Methods

第十四章: 线性规划之内点法

上世纪八十年代, 人们发现 许多大型的线性规划问题可用非线性规划和非线性方程组的构造和算法高效地求解. 这些方法的一个特征是,
$$它们均需要所有的迭代点\underline{\text{严格满足}}问题中的不等式约束(后面也发展了不可行方法),$$因此这些方法也被称作 内点法(interior-point methods). 在九十年代早期以前, 内点法的一个子类—— 原始-对偶算法(primal-dual methods)——脱颖而出, 成为最高效的实用算法. 其 在大型问题上的表现可与单纯形法相媲美. 本章我们就聚焦原始-对偶算法.

内点法出现来源于对理论性质优于单纯形法算法的需求. 正如我们在第十三章谈到的, 单纯形法在特定的病态问题上可能会十分低效. 粗略地讲, 使用单纯形法求解一个线性规划的时间可能会达问题大小的指数级别. 这里问题的大小由未知量的数量和所需的储存量所度量. 尽管对于大多数实际问题, 单纯形法的效果都要远远好于这个极端, 但它这一最坏情形复杂度却推动了对具有更佳表现算法的发展. 第一类这样的算法是Khachiyan提出的椭圆方法(ellipsoid method), 它能保证在多项式时间内找到问题的解. 但不幸的是, 此法在所有的问题上均接近以其最差情形表现, 因此在实际应用中还无法与单纯形法竞争.
Karmarkar在1984年提出的投影算法projective algorithm同样具有多项式复杂度, 但它因良好的实际表现而备受关注. 尽管最初关于它在大规模线性规划上具有出色表现的断言从未被完全证实, 但这的确推动了许多新方法的发展. 这些新方法均与Karmarkar最初的算法和之后要讨论的对数障碍函数法log-barrier approach相关, 不过事实上它们中的许多都可以由早期的方法独立地推导出来.

每种内点法与单纯形法都有一样的不同点:

• 内点法的每步迭代计算代价较大, 但这些迭代步可以很快地收敛到解; 单纯形法通常需要大量经济却较低效的迭代.
• 从几何上看, 单纯形法的搜索路径始终位于可行多面体的边界, 而内点法从内部或外部趋近于位于边界的解.

本章的内容大致如下:

1. 我们将概述原始-对偶内点法背后的一些基本想法, 包括它们与牛顿法和同伦(连续)方法的关系, 并介绍中心路径的概念.
2. 我们将粗略地介绍这类方法中较为重要的成员, 并对一种称为长步路径跟踪算法进行全面的收敛分析.
3. 我们会详细地介绍Mehrotra提出的一种实用的预估-校正算法. 它是当前许多软件的基础.

1. 原始-对偶算法

1.1 概述

我们考虑标准形式的线性规划问题:$$\min c^{T}x,\mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}Ax=b,x\ge 0,$$其中$c,x\in {\mathbb{R}}^n,b\in {\mathbb{R}}^m$, $A$为$m\times n$的行满秩矩阵. 这一问题的对偶问题为$$\max b^T\lambda ,\mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}{A}^T\lambda +s=c,s\ge 0,$$其中$\lambda \in {\mathbb{R}}^m,s\in {\mathbb{R}}^n$. 如第十三章讨论, 原始和对偶问题的解由以下KKT条件所刻画:$$\begin{array}{rl}{A}^T\lambda +s& =c,\\ Ax& =b,\\ x_i{s}_i& =0,i=1,2,\dots ,n,\\ (x,s)& \ge 0.\end{array}$$原始-对偶算法通过将Newton法的变体形式应用于上述前三个等式, 并修正搜索方向和步长以满足严格不等式条件, 从而求满足KKT条件的解$(x^{\ast },{\lambda }^{\ast },s^{\ast })$. 注意前三个等式条件要么是线性的, 要么稍微有点非线性, 因此单就求解它们而言并不困难. 但若加入对不等式条件的考虑, 情况就要复杂得多. 这也给内点法的设计和分析带来了困难.

为推导原始-对偶内点法, 我们以另一种稍微不同的形式重述最优性条件: $$\begin{array}{rl}F(x,\lambda ,s)=\left[\begin{array}{c}{A}^T\lambda +s-c\\ Ax-b\\ XSe\end{array}\right]& =0,\\ (x,s)& \ge 0,\end{array}$$其中$$\begin{gathered} F:{\mathbb{R}}^{2n+m}\to {\mathbb{R}}^{2n+m}, \\ X=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(x_1,x_2,\dots ,x_n),S=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(s_1,s_2,\dots ,s_n), \\ e=(1,1,\dots ,1{)}^T. \end{gathered}$$原始-对偶法将产生严格满足界约束的迭代点$(x^k,{\lambda }^k,s^k):x^k>0,s^k>0$. 这一性质是"内点法"这一称谓的来源. 通过对界约束的要求, 算法就能规避假解, 即那些满足$F(x,\lambda ,s)=0$但却不满足$(x,s)\ge 0$的解. 假解有很多, 并且并不能提供关于问题解的有效信息, 因此将它们排除出搜索区域是很有意义的.

如优化中大多数迭代算法一样, 原始-对偶内点法有两个基本要素:

1. 确定迭代步的程序. 它根植于求解非线性方程组的牛顿法. 它在当前点附近构造$F$的线性模型, 通过求解以下线性方程组得到搜索方向$(\Delta x,\Delta \lambda ,\Delta s)$:$$J(x,\lambda ,s)\left[\begin{array}{c}\Delta x\\ \Delta \lambda \\ \Delta s\end{array}\right]=-F(x,\lambda ,s),$$其中$J$为$F$的Jacobi矩阵. 关于非线性方程组的牛顿法的讨论可见第十一章. 记$$r_b=Ax-b,r_c=A^T\lambda +s-c,$$则可将牛顿方程写成$$\left[\begin{array}{ccc}O& A^T& I\\ A& O& O\\ S& O& X\end{array}\right]\left[\begin{array}{c}\Delta x\\ \Delta \lambda \\ \Delta s\end{array}\right]=\left[\begin{array}{c}-r_c\\ -r_b\\ -XSe\end{array}\right].$$
2. 在搜索空间中对每个点好坏与否的度量. 这里度量需要包含一个很重要的部分, 即在$x>0,s>0$时$x_i{s}_i,i=1,2,\dots ,n$的均值. 这个量被称为对偶度量(duality measure):$$\mu =\frac{1}{n}\sum _{i=1}^n{x}_i{s}_i=\frac{x^{T}s}{n}.$$

一般地, 沿着纯牛顿方向走一满步(full step)是会违反界约束$(x,s)\ge 0$的, 因此我们需要沿着牛顿方向做一次线搜索并定义新的迭代点为$$(x,\lambda ,s)+\alpha (\Delta x,\Delta \lambda .\Delta s),$$其中$\alpha \in (0,1]$为某个符合线搜索条件的步长. 通常我们在打破条件$(x,s)>0$前只能沿着牛顿方向走很小的一步($\alpha \ll 1$). 纯粹的牛顿方向(有时也被称为仿射变换方向(affine scaling direction))往往并不能给我们带来迭代点上多大的改进.

正是因此, 大多原始-对偶算法使用更加"温和"的牛顿方向, 即算法并不直接求解牛顿方程, 而是抱着降低对偶度量(不是到0)的目的求解一类修正方程. 特别地, 我们想得到牛顿步使得下一步的$x_i{s}_i$趋近于$x_i{s}_i=\sigma \mu$, 其中$\mu$为当前的对偶度量, $\sigma \in [0,1]$为我们希望在此步缩减对偶度量的比例. 于是修正的牛顿方程为$$\left[\begin{array}{ccc}O& A^T& I\\ A& O& O\\ S& O& X\end{array}\right]\left[\begin{array}{c}\Delta x\\ \Delta \lambda \\ \Delta s\end{array}\right]=\left[\begin{array}{c}-r_c\\ -r_b\\ -XSe+\sigma \mu e\end{array}\right].$$事实上, 从第三行可知$$S\Delta x+X\Delta s=-XSe+\sigma \mu e.$$因此$$\begin{array}{rl}(x+\Delta x{)}^T(s+\Delta s)& =x^{T}s+x^T\Delta s+s^T\Delta x+\Delta x^T\Delta s\\ & =n\mu +e^{T}X\Delta s+e^{T}S\Delta x+\Delta x^T\Delta s\\ & =n\mu +e^T(-XSe+\sigma \mu e)+\Delta x^T\Delta s\\ & =n\sigma \mu +\Delta x^T\Delta s.\end{array}$$当$\Delta x,\Delta s$为小量时, 新迭代点的对偶度量就近似等于$\sigma \mu$, 于是每个$x_i{s}_i\approx \sigma \mu$. 我们称$\sigma$为中心化参数(centering parameter), 这一称谓我们后面会道出其来源. 当$\sigma >0$, 算法往往可以在不打破$(x,s)\ge 0$的前提下走更长的$\alpha$. 至今, 我们已介绍了路径跟踪path-following原始-对偶内点法的大部内容. 此类算法的一般框架如下.

框架1 (Primal-Dual Path-Following)
Given $(x^0,{\lambda }^0,s^0)$ with $(x^0,s^0)>0$;
for $k=0,1,2,\dots$
\quad\quad Choose ${\sigma }_k\in [0,1]$ and solve$$\left[\begin{array}{ccc}O& A^T& I\\ A& O& O\\ S^k& O& X^k\end{array}\right]\left[\begin{array}{c}\Delta x^k\\ \Delta {\lambda }^k\\ \Delta s^k\end{array}\right]=\left[\begin{array}{c}-r_c^k\\ -r_b^k\\ -X^k{S}^{k}e+{\sigma }_k{\mu }_{k}e\end{array}\right],$$ \quad\quad\quad\quad where ${\mu }_k=(x^k{)}^T{s}^k/n$;
\quad\quad Set$$(x^{k+1},{\lambda }^{k+1},s^{k+1})=(x^k,{\lambda }^k,s^k)+{\alpha }_k(\Delta x^k,\Delta {\lambda }^k,\Delta s^k),$$ \quad\quad\quad\quad choosing ${\alpha }_k$ so that $(x^{k+1},s^{k+1})>0$.
end (for)

中心化参数${\sigma }_k$和步长${\alpha }_k$的选取对于算法的表现至关重要. 直接或间接控制这些参数的技术将得到具有不同性质的方法.

尽管实施内点法的软件通常并不从一可行点$(x^0,{\lambda }^0,s^0)$(这里"可行"相应于$r_b̸=0$或$r_c̸=0$)出发, 但大多理论和算法中均假设这些条件是满足的. 在本节剩下的部分, 我们讨论可行方法, 使用基础的数学工具和概念, 进行全面的收敛性分析. 不可行方法的分析套路是类似的, 但其细节要更加复杂, 这里就不多说明. 不过在第2节, 我们会介绍一个不需要可行初始点的完整实用算法.
为开启我们的讨论和分析, 我们首先引入中心路径central path和它邻域的概念.

1.2 中心路径

我们定义原始-对偶可行集$\mathcal{F}$和严格可行集${\mathcal{F}}^o$为 F = { ( x , λ , s ) ∣ A x = b , A T λ + s = c , ( x , s ) ≥ 0 } , F o = { ( x , λ , s ) ∣ A x = b , A T λ + s = c , ( x , s ) > 0 } . \begin{aligned}\mathcal{F}&=\{(x,\lambda,s)\mid Ax=b,A^T\lambda+s=c,(x,s)\ge0\},\\\mathcal{F}^o&=\{(x,\lambda,s)\mid Ax=b,A^T\lambda+s=c,(x,s)>0\}.\end{aligned} FFo​={(x,λ,s)∣Ax=b,ATλ+s=c,(x,s)≥0},={(x,λ,s)∣Ax=b,ATλ+s=c,(x,s)>0}.​中心路径$\mathcal{C}$为一条由严格可行点构成的弧段, 它在原始-对偶算法中扮演着重要的角色. 它的参数为$\tau >0$, $\mathcal{C}$上的每一点$(x_{\tau },{\lambda }_{\tau },s_{\tau })$满足以下方程:$$\begin{array}{rl}{A}^T\lambda +s& =c,\\ Ax& =b,\\ x_i{s}_i& =\tau ,i=1,2,\dots ,n,\\ (x,s)& >0.\end{array}$$上述方程与KKT条件仅在互补松弛条件的右端项上不同. 这里我们要求每个乘积$x_i{s}_i$取相同的正值. 于是中心路径就有参数化表示 C = { ( x τ , λ τ , s τ ) ∣ τ > 0 } . \mathcal{C}=\{(x_{\tau},\lambda_{\tau},s_{\tau})\mid\tau>0\}. C={(xτ​,λτ​,sτ​)∣τ>0}.可以证明对每个$\tau >0$, $(x_{\tau },{\lambda }_{\tau },s_{\tau })$是唯一确定的当且仅当${\mathcal{F}}^o$非空.

定义中心路径的方式还有:

1. 原始问题的对数障碍logarithmic-barrier重构问题的最优性条件. 通过对非负约束引入对数障碍(log-barrier)项(障碍参数$\tau >0$), 我们得到$$\min c^{T}x-\tau \sum _{i=1}^n\ln x_i,\mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}Ax=b.$$此问题的KKT条件为$$c_i-\frac{\tau }{x_i}-A_{\cdot i}^T\lambda =0,i=1,2,\dots ,n;Ax=b,$$其中$\lambda$为对应于等式约束的Lagrange乘子. 由于目标函数严格凸, 因此KKT条件充分必要. 若定义$s_i=\tau /x_i,i=1,2,\dots ,n$, 则就回到了中心路径满足的方程.

2. 利用映射$F:$$$F(x_{\tau },{\lambda }_{\tau },s_{\tau })=\left[\begin{array}{c}0\\ 0\\ \tau e\end{array}\right],(x_{\tau },s_{\tau })>0.$$

中心路径所满足的参数方程在$\tau \to 0$时将逐渐逼近原始问题的KKT条件. 倘若$\mathcal{C}$随着$\tau \downarrow 0$收敛, 则必定收敛到线性规划的原始-对偶解. 因此中心路径就是一条保持$x,s$严格正性且同步降低每个乘积$x_i{s}_i,i=1,2,\dots ,n$的轨迹.

大多原始-对偶算法都采用指向$\mathcal{C}$的牛顿步(即$\tau >0$), 而不是$F$的纯牛顿步. 指向$\mathcal{C}$的牛顿步会带来的好处是, 它们往往指向非负约束构作的非负象限的内部, 因此沿着下一步搜索方向搜索比沿着纯牛顿步(仿射变换)能走得更远.

当$(x,\lambda ,s)\in \mathcal{F}$时, 我们有$r_b=0,r_c=0$, 因此搜索方向满足$$\left[\begin{array}{ccc}O& A^T& I\\ A& O& O\\ S& O& X\end{array}\right]\left[\begin{array}{c}\Delta x\\ \Delta \lambda \\ \Delta s\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ -XSe+\sigma \mu e\end{array}\right],$$其中$\mu$为对偶度量, $\sigma \in [0,1]$为中心化参数. 特别地,

1. 当$\sigma =1$时, 上述方程的解称为中心化方向(centering direction), 即指向$\mathcal{C}$上的一点$(x_{\mu },{\lambda }_{\mu },s_{\mu }):x_i{s}_i=\mu$. 中心化方向通常严格指向非负象限内部, 并且往往并不会给对偶度量$\mu$带来多大的改进. 不过, 由于离得$\mathcal{C}$更近了, 我们就可期待下一次迭代能带来的改进.
2. 当$\sigma =0$时, 这就是标准的牛顿步(仿射变换步).

这就是为什么$\sigma$被称作中心化参数的缘由. 许多算法使用开区间$(0,1)$中的$\sigma$以权衡减小$\mu$和改善中心性的双重目标.

1.3 中心路径的邻域和路径跟踪算法

路径跟随(path-following)算法显式地将迭代点限制在中心路径$\mathcal{C}$的一个邻域内, 并随着$\mathcal{C}$收敛到线性规划的解. 通过防止迭代点过于接近非负象限的边界, 它们保证每一步迭代均是非平凡的. 进一步地, 通过限制对偶度量${\mu }_k\to 0(k\to \infty )$, 我们可保证迭代点$(x^k,{\lambda }^k,s^k)$越来越趋近于满足KKT条件.

1.3.1 中心路径的邻域

$\mathcal{C}$的两种邻域为 N 2 ( θ ) = { ( x , λ , s ) ∈ F o ∣ ∥ X S e − μ e ∥ 2 ≤ θ μ } , θ ∈ [ 0 , 1 ) , \mathcal{N}_2(\theta)=\{(x,\lambda,s)\in\mathcal{F}^o\mid\Vert XSe-\mu e\Vert_2\le\theta\mu\},\quad \theta\in[0,1), N2​(θ)={(x,λ,s)∈Fo∣∥XSe−μe∥2​≤θμ},θ∈[0,1),以及 N − ∞ ( γ ) = { ( x , λ , s ) ∈ F o ∣ x i s i ≥ γ μ , i = 1 , 2 , … , n } , γ ∈ ( 0 , 1 ] . \mathcal{N}_{-\infty}(\gamma)=\{(x,\lambda,s)\in\mathcal{F}^o\mid x_is_i\ge\gamma\mu,i=1,2,\ldots,n\},\quad \gamma\in(0,1]. N−∞​(γ)={(x,λ,s)∈Fo∣xi​si​≥γμ,i=1,2,…,n},γ∈(0,1].
这两个邻域参数的常用取值为$\theta =0.5,\gamma =10^{-3}$. 若一点位于${\mathcal{N}}_{-\infty }(\gamma )$, 则每个乘积$x_i{s}_i$至少为它们均值$\mu$的$\gamma$倍. 这个条件并不高, 并且我们可以通过令$\gamma \to 0$, 扩充${\mathcal{N}}_{-\infty }(\gamma )$至整个可行集$\mathcal{F}$. 相比之下, ${\mathcal{N}}_2(\theta )$的限制就大一些, 毕竟不论$\theta$不论有多么接近1, 总可能有${\mathcal{F}}^o$中的点不在${\mathcal{N}}_2(\theta )$中.

通过限制所有的迭代点于邻域中, 路径跟踪算法就差不多是以同样的速率减小所有的乘积$x_i{s}_i$. 下图就展示了将中心路径$\mathcal{C}$投影到原始变量空间的情形, 同时还可见$\mathcal{C}$的邻域$\mathcal{N}$.

1.3.2 路径跟踪算法

路径跟踪算法与处理一般非线性方程组的同伦方法很相似, 后者也定义了一条通向解的路径. 二者不同在于:

• 传统的同伦算法待在路径的一个管状邻域内, 增大参数, 跟随路径直至解.
• 对于原始-对偶算法, 中心路径的邻域则是角状的而不是管状的: 它随$\mu$的增大而变得宽厚松弛, 随$\mu$的减小而变得窄小萎缩(这是因为要保证$(x,s)>0$).

下面的算法是框架1的特殊情形, 被称作是长步路径跟踪算法(long-step path-following algorithm). 此算法因使用邻域${\mathcal{N}}_{-\infty }(\gamma )$而获得了较快的收敛速度($\gamma \to 0$). 它还需要两个参数——${\sigma }_{\min },{\sigma }_{\max }$——分别是中心化参数${\sigma }_k$的下界和上界. 通常, 搜索方向是通过求解修正牛顿方程得到的, 并且我们将在不逾越邻域边界的前提下, 选取尽可能大的步长${\alpha }_k$.
此后的分析中, 我们使用记号:$$\begin{array}{rl}(x^k(\alpha ),{\lambda }^k(\alpha ),s^k(\alpha ))& \stackrel{def}{=}(x^k,{\lambda }^k,s^k)+\alpha (\Delta x^k,\Delta {\lambda }^k,\Delta s^k),\\ {\mu }_k(\alpha )& \stackrel{def}{=}{x}^k(\alpha {)}^T{s}^k(\alpha )/n.\end{array}$$

算法2 (Long-Step Path-Following)
Given $\gamma ,{\sigma }_{\min },{\sigma }_{\max }$ with $\gamma \in (0,1),0<{\sigma }_{\min }\le {\sigma }_{\max }<1$, and $(x^0,{\lambda }^0,s^0)\in {\mathcal{N}}_{-\infty }(\gamma )$;
for $k=0,1,2,\dots$
\quad\quad Choose ${\sigma }_k\in [{\sigma }_{\min },{\sigma }_{\max }]$;
\quad\quad Solve the modified Newton equations to obtain $(\Delta x^k,\Delta {\lambda }^k,\Delta s^k)$;
\quad\quad Choose ${\alpha }_k$ as the largest value of $\alpha$ in $[0,1]$ such that $$(x^k(\alpha ),{\lambda }^k(\alpha ),s^k(\alpha ))\in {\mathcal{N}}_{-\infty }(\gamma );$$
\quad\quad Set $(x^{k+1},{\lambda }^{k+1},s^{k+1})=(x^k({\alpha }_k),{\lambda }^k({\alpha }_k),s^k({\alpha }_k))$;
end (for)

注意这里选取步长时不再像之前那样需要用满足一定条件(如Wolfe条件)的线搜索. 一方面我们缺少选取步长的价值函数, 另一方面中心路径已经保证了收敛.

算法的表现可见下图$n=2$的例子.

上图中横纵坐标分别表示乘积$x_1{s}_1,x_2{s}_2$, 因此中心路径$\mathcal{C}$为从原点出发与横轴夹$45^{\circ }$角的射线. (这里的原点若满足可行性条件则就是一个原始-对偶解.) 搜索方向则以曲线呈现.

上图表示, 中心化参数的下界${\sigma }_{\min }$使得每个搜索方向都先远离${\mathcal{N}}_{-\infty }(\gamma )$的边界进入邻域的相对内部. 即, 沿搜索方向的一小步可以改善中心性. 选取较大的步长$\alpha$则又会逾越邻域的边界, 这是因为当$\alpha$增加时, 使用线性系统近似非线性系统的误差将会增大. 不过, 我们仍然可以保证在到达${\mathcal{N}}_{-\infty }$的边界之前我们是可以走一小步的.

下面我们将分析算法2. 通过谨慎的选取${\sigma }_k$, 算法2在实际使用时是十分高效的. 只需经过一些细节上的完善, 它就会成为一个真正强大的算法. 这一点我们将在下一节讨论之.

我们将要证明: 给定容忍限$ϵ>0$, 算法需要$O(n|\log ϵ|)$步迭代将对偶度量减小为原来的$ϵ$倍, 即寻得一点$(x^k,{\lambda }^k,s^k):{\mu }_k\le ϵ{\mu }_0$. 对于充分小的$ϵ$, $(x^k,{\lambda }^k,s^k)$就越接近于原始-对偶解. 而迭代步数的估计$O(n|\log ϵ|)$的估计仅仅是最差情形估计, 在实际问题中, $n$的增加仅会轻微改变所需的迭代数. 这就像单纯形法, 对于Klee-Minty反例, 它需要$2^n$次迭代才能得到解, 而一般它只需$O(m)$次迭代, 这里$m$为约束矩阵$A$的行数.

下面我们先证明一个技术性的引理(引理1)并使用它推导乘积$\Delta x_i\Delta s_i$构成向量的界(引理2). 接着定理3将寻得步长${\alpha }_k$的下界和对应的在第$k$次迭代上$\mu$的下降估计. 最终定理4证明对于$ϵ\in (0,1)$, 若要${\mu }_k<ϵ$至多需要经过$O(n|\log ϵ|)$次迭代.

引理1 $u,v\in {\mathbb{R}}^n:u^{T}v\ge 0$, 则$$\Vert UVe{\Vert }_2\le 2^{-3/2}\Vert u+v{\Vert }_2^2,$$其中$$U=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(u_1,u_2,\dots ,u_n),V=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(v_1,v_2,\dots ,v_n).$$
证明: 首先对于任意两个标量$\alpha ,\beta :\alpha \beta \ge 0$, 我们有算术-几何平均值不等式$$\sqrt{|\alpha \beta |}\le \frac{1}{2}|\alpha +\beta |.$$由于$u^{T}v\ge 0$, 因此$$0\le u^{T}v=\sum _{u_i{v}_i\ge 0}{u}_i{v}_i+\sum _{u_i{v}_i<0}{u}_i{v}_i=\sum _{i\in \mathcal{P}}|u_i{v}_i|-\sum _{i\in \mathcal{M}}|u_i{v}_i|,$$这里我们将指标集 { 1 , 2 , … , n } \{1,2,\ldots,n\} {1,2,…,n}划分为 P = { i ∣ u i v i ≥ 0 } , M = { i ∣ u i v i &lt; 0 } . \mathcal{P}=\{i\mid u_iv_i\ge0\},\quad \mathcal{M}=\{i\mid u_iv_i&lt;0\}. P={i∣ui​vi​≥0},M={i∣ui​vi​<0}.于是$$\begin{array}{rl}\Vert UVe\Vert & ={\left(\Vert [u_i{v}_i{]}_{i\in \mathcal{P}}{\Vert }^2+\Vert [u_i{v}_i{]}_{i\in \mathcal{M}}{\Vert }^2\right)}^{1/2}\\ & \le {\left(\Vert [u_i{v}_i{]}_{i\in \mathcal{P}}{\Vert }_1^2+\Vert [u_i{v}_i{]}_{i\in \mathcal{M}}{\Vert }_1^2\right)}^{1/2}(\because \Vert \cdot {\Vert }_2\le \Vert \cdot {\Vert }_1)\\ & \le {\left(2\Vert [u_i{v}_i{]}_{i\in \mathcal{P}}{\Vert }_1^2\right)}^{1/2}\\ & \le \sqrt{2}{\Vert {\left[\frac{1}{4}(u_i+v_i{)}^2\right]}_{i\in \mathcal{P}}\Vert }_1\\ & =2^{-3/2}\sum _{i\in \mathcal{P}}(u_i+v_i{)}^2\\ & \le 2^{-3/2}\sum _{i=1}^n(u_i+v_i{)}^2\\ & =2^{-3/2}\Vert u+v{\Vert }^2,\end{array}$$证毕.

下面的结论和证明中, 我们省去下标$k$并类似地定义对角阵$\Delta X,\Delta S$为$$\Delta X=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(\Delta x_1,\Delta x_2,\dots ,\Delta x_n),\Delta S=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(\Delta s_1,\Delta s_2,\dots ,\Delta s_n).$$

引理2 对$(x,\lambda ,s)\in {\mathcal{N}}_{-\infty }(\gamma )$, 有$$\Vert \Delta X\Delta Se\Vert \le 2^{-3/2}(1+1/\gamma )n\mu .$$
证明: 显然应有$\Delta x^T\Delta s=0$. 在修正牛顿方程的最后一行左乘$(XS{)}^{-1/2}$, 利用对角阵相乘可交换, 可得$$D^{-1}\Delta x+D\Delta s=(XS{)}^{-1/2}(-XSe+\sigma \mu e).$$因$(D^{-1}\Delta x{)}^T(D\Delta s)=\Delta x^T\Delta s=0$, 因此我们可将引理1应用于$u=D^{-1}\Delta x,v=D\Delta s$得到$$\begin{array}{rl}\Vert \Delta X\Delta Se\Vert & =\Vert (D^{-1}\Delta X)(D\Delta S)e\Vert \\ & \le 2^{-3/2}\Vert D^{-1}\Delta x+D\Delta s{\Vert }^2(由引理1)\\ & =2^{-3/2}\Vert (XS{)}^{-1/2}(-XSe+\sigma \mu e){\Vert }^2\\ & =2^{-3/2}\left[x^{T}s-2\sigma \mu e^{T}e+{\sigma }^2{\mu }^2\sum _{i=1}^n\frac{1}{x_i{s}_i}\right]\\ & \le 2^{-3/2}\left[x^{T}s-2\sigma \mu e^{T}e+{\sigma }^2{\mu }^2\frac{n}{\gamma \mu }\right]\\ & \le 2^{-3/2}\left[1-2\sigma +\frac{{\sigma }^2}{\gamma }\right]n\mu \\ & \le 2^{-3/2}(1+1/\gamma )n\mu .\end{array}$$

定理3 对于算法2, 给定参数$\gamma ,{\sigma }_{\min },{\sigma }_{\max }$, 存在不依赖$n$的常数$\delta$使得$${\mu }_{k+1}\le \left(1-\frac{\delta }{n}\right){\mu }_k,\forall k\ge 0.$$
证明: 我们先证明$$(x^k(\alpha ),{\lambda }^k(\alpha ),s^k(\alpha ))\in {\mathcal{N}}_{-\infty }(\gamma ),\forall \alpha \in \left[0,2^{3/2}\gamma \frac{1-\gamma }{1+\gamma }\frac{{\sigma }_k}{n}\right].$$于是立刻就可得步长${\alpha }_k$应至少是此区间的上界, 即$${\alpha }_k\ge 2^{3/2}\gamma \frac{1-\gamma }{1+\gamma }\frac{{\sigma }_k}{n}.$$对$i=1,2,\dots ,n$, 从引理2可得$$|\Delta x_i^k\Delta s_i^k|\le \Vert \Delta X^k\Delta S^{k}e{\Vert }_2\le 2^{-3/2}(1+1/\gamma )n{\mu }_k.$$由$x_i^k{s}_i^k\ge \gamma {\mu }_k$得到$$\begin{array}{rl}{x}_i^k(\alpha )s_i^k(\alpha )& =(x_i^k+\alpha \Delta x_i^k)(s_i^k+\alpha \Delta s_i^k)\\ & =x_i^k{s}_i^k+\alpha (x_i^k\Delta s_i^k+s_i^k\Delta x_i^k)+{\alpha }^2\Delta x_i^k\Delta s_i^k\\ & \ge x_i^k{s}_i^k(1-\alpha )+\alpha {\sigma }_k{\mu }_k-{\alpha }^2|\Delta x_i^k\Delta s_i^k|\\ & \ge \gamma (1-\alpha ){\mu }_k+\alpha {\sigma }_k{\mu }_k-{\alpha }^2{2}^{-3/2}(1+1/\gamma )n{\mu }_k.\end{array}$$将$S^k\Delta x^k+X^k\Delta s^k=-X^k{S}^{k}e+{\sigma }_k{\mu }_{k}e$的$n$个分量相加并利用${\mu }_k,{\mu }_k(\alpha )$的定义, 可得$${\mu }_k(\alpha )=(1-\alpha (1-{\sigma }_k)){\mu }_k.$$由上面两个式子, 若再有$$\gamma (1-\alpha ){\mu }_k+\alpha {\sigma }_k{\mu }_k-{\alpha }^2{2}^{-3/2}(1+1/\gamma )n{\mu }_k\ge \gamma (1-\alpha +\alpha {\sigma }_k){\mu }_k,$$则就有近似条件proximity condition成立:$$x_i^k(\alpha )s_i^k(\alpha )\ge \gamma {\mu }_k(\alpha ).$$整理前面的条件可得$$\begin{gathered} \alpha {\sigma }_k{\mu }_k(1-\gamma )\ge {\alpha }^2{2}^{-3/2}n{\mu }_k(1+1/\gamma ) \\ \Rightarrow \alpha \le \frac{2^{3/2}}{n}{\sigma }_k\gamma \frac{1-\gamma }{1+\gamma }. \end{gathered}$$
下面我们估计在第$k$步对偶度量$\mu$的下降. $$\begin{array}{rl}{\mu }_{k+1}& =x^k({\alpha }_k{)}^T{s}^k({\alpha }_k)/n\\ & =\left[(x^k{)}^T{s}^k+\alpha \left((x^k{)}^T\Delta s^k+(s^k{)}^T\Delta x^k\right)+{\alpha }_k^2(\Delta x^k{)}^T\Delta s^k\right]/n\\ & ={\mu }_k+{\alpha }_k\left(-(x^k{)}^T{s}^k/n+{\sigma }_k{\mu }_k\right)\\ & =(1-{\alpha }_k(1-{\sigma }_k)){\mu }_k\\ & \le \left(1-\frac{2^{3/2}}{n}\gamma \frac{1-\gamma }{1+\gamma }{\sigma }_k(1-{\sigma }_k)\right){\mu }_k.\end{array}$$因$\sigma (1-\sigma )$为凹二次函数, 因此其在任何区间上的最小值必在端点处取到. 因此我们有 σ k ( 1 − σ k ) ≥ min ⁡ { σ min ⁡ ( 1 − σ min ⁡ ) , σ max ⁡ ( 1 − σ max ⁡ ) } , ∀ σ k ∈ [ σ min ⁡ , σ max ⁡ ] . \sigma_k(1-\sigma_k)\ge\min\{\sigma_{\min}(1-\sigma_{\min}),\sigma_{\max}(1-\sigma_{\max})\},\quad\forall\sigma_k\in[\sigma_{\min},\sigma_{\max}]. σk​(1−σk​)≥min{σmin​(1−σmin​),σmax​(1−σmax​)},∀σk​∈[σmin​,σmax​].代入并置 δ = 2 3 / 2 γ 1 − γ 1 + γ min ⁡ { σ min ⁡ ( 1 − σ min ⁡ ) , σ max ⁡ ( 1 − σ max ⁡ ) } . \delta=2^{3/2}\gamma\frac{1-\gamma}{1+\gamma}\min\{\sigma_{\min}(1-\sigma_{\min}),\sigma_{\max}(1-\sigma_{\max})\}. δ=23/2γ1+γ1−γ​min{σmin​(1−σmin​),σmax​(1−σmax​)}.证毕.

最后我们证明对偶度量$\mu$下降为原来的$ϵ$倍至多需要$O(n\log 1/ϵ)$次迭代.

定理4 给定$ϵ\in (0,1),\gamma \in (0,1)$, 假定算法2的初始点满足$(x^0,{\lambda }^0,s^0)\in {\mathcal{N}}_{-\infty }(\gamma )$. 于是存在$K=O(n\log 1/ϵ)$使得$${\mu }_k\le ϵ{\mu }_0,\forall k\ge K.$$
证明: 对定理3的结论取对数可得$$\log {\mu }_{k+1}\le \log \left(1-\frac{\delta }{n}\right)+\log {\mu }_k.$$重复迭代, 有$$\log {\mu }_k\le k\log \left(1-\frac{\delta }{n}\right)+\log {\mu }_0.$$从不等式估计$$\log (1+\beta )\le \beta ,\forall \beta >-1$$可知$$\log ({\mu }_k/{\mu }_0)\le k\left(-\frac{\delta }{n}\right).$$因此${\mu }_k/{\mu }_0\le ϵ$成立当$$k\left(-\frac{\delta }{n}\right)\le \log ϵ.$$自然这还需要$\delta /n<1$. 此不等式对$$k\ge K\stackrel{def}{=}\frac{n}{\delta }\log \frac{1}{ϵ}=\frac{n}{\delta }|\log ϵ|$$成立, 证毕.

2. 实用原始-对偶算法

实用内点法的实施继承了之前一小节的内容: 算法过程保证$x^k,s^k$的严格正性, 且每一步都涉及中心化的类牛顿步. 不过, 大多数算法实施不论是初始点还是在迭代过程中都是不可行的. 本节, 我们忽略此类不可行算法的理论分析, 而仅讨论在实用层面上具有重大效用的强化方法, 最终归总于算法3.

2.1 校正步和中心化步

实用算法的一个重要特征是, 它们使用校正步corrector steps来弥补由线性化非线性项$x_i{s}_i=0$所带来的误差. 考虑仿射变换方向$(\Delta x,\Delta \lambda ,\Delta s):$$$\left[\begin{array}{ccc}O& A^T& I\\ A& O& O\\ S& O& X\end{array}\right]\left[\begin{array}{c}\Delta x^{\mathrm{a}\mathrm{f}\mathrm{f}}\\ \Delta {\lambda }^{\mathrm{a}\mathrm{f}\mathrm{f}}\\ \Delta s^{\mathrm{a}\mathrm{f}\mathrm{f}}\end{array}\right]=\left[\begin{array}{c}-r_c\\ -r_b\\ -XSe\end{array}\right].$$若沿着此方向走一满步, 则有$$\begin{array}{rl}& (x_i+\Delta x_i^{\mathrm{a}\mathrm{f}\mathrm{f}})(s_i+\Delta s_i^{\mathrm{a}\mathrm{f}\mathrm{f}})\\ & =x_i{s}_i+x_i\Delta s_i^{\mathrm{a}\mathrm{f}\mathrm{f}}+s_i\Delta x_i^{\mathrm{a}\mathrm{f}\mathrm{f}}+\Delta x_i^{\mathrm{a}\mathrm{f}\mathrm{f}}\Delta s_i^{\mathrm{a}\mathrm{f}\mathrm{f}}=\Delta x_i^{\mathrm{a}\mathrm{f}\mathrm{f}}\Delta s_i^{\mathrm{a}\mathrm{f}\mathrm{f}}.\end{array}$$这就是说, 更新后的$x_i{s}_i$不是0而是$\Delta x_i^{\mathrm{a}\mathrm{f}\mathrm{f}}\Delta s_i^{\mathrm{a}\mathrm{f}\mathrm{f}}$. 我们可以通过求解以下系统得到$(\Delta x^{\mathrm{c}\mathrm{o}\mathrm{r}},\Delta {\lambda }^{\mathrm{c}\mathrm{o}\mathrm{r}},\Delta s^{\mathrm{c}\mathrm{o}\mathrm{r}})$以校正这一偏差:$$\left[\begin{array}{ccc}O& A^T& I\\ A& O& O\\ S& O& X\end{array}\right]\left[\begin{array}{c}\Delta x^{\mathrm{c}\mathrm{o}\mathrm{r}}\\ \Delta {\lambda }^{\mathrm{c}\mathrm{o}\mathrm{r}}\\ \Delta s^{\mathrm{c}\mathrm{o}\mathrm{r}}\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ -\Delta X^{\mathrm{a}\mathrm{f}\mathrm{f}}\Delta S^{\mathrm{a}\mathrm{f}\mathrm{f}}e\end{array}\right].$$这样, $$\begin{array}{rl}& (x_i+\Delta x_i^{\mathrm{a}\mathrm{f}\mathrm{f}}+\Delta x_i^{\mathrm{c}\mathrm{o}\mathrm{r}})(s_i+\Delta s_i^{\mathrm{a}\mathrm{f}\mathrm{f}}+\Delta s_i^{\mathrm{c}\mathrm{o}\mathrm{r}})\\ & =\Delta x_i^{\mathrm{a}\mathrm{f}\mathrm{f}}\Delta s_i^{\mathrm{a}\mathrm{f}\mathrm{f}}+(x_i+\Delta x_i^{\mathrm{a}\mathrm{f}\mathrm{f}})\Delta s_i^{\mathrm{c}\mathrm{o}\mathrm{r}}+\Delta x_i^{\mathrm{c}\mathrm{o}\mathrm{r}}(s_i+\Delta s_i^{\mathrm{a}\mathrm{f}\mathrm{f}})+\Delta x_i^{\mathrm{c}\mathrm{o}\mathrm{r}}\Delta s_i^{\mathrm{c}\mathrm{o}\mathrm{r}}\\ & =\Delta x_i^{\mathrm{a}\mathrm{f}\mathrm{f}}\Delta s_i^{\mathrm{a}\mathrm{f}\mathrm{f}}-\Delta x_i^{\mathrm{a}\mathrm{f}\mathrm{f}}\Delta s_i^{\mathrm{a}\mathrm{f}\mathrm{f}}+\Delta x_i^{\mathrm{a}\mathrm{f}\mathrm{f}}\Delta s_i^{\mathrm{c}\mathrm{o}\mathrm{r}}+\Delta x_i^{\mathrm{c}\mathrm{o}\mathrm{r}}\Delta s_i^{\mathrm{a}\mathrm{f}\mathrm{f}}+\Delta x_i^{\mathrm{c}\mathrm{o}\mathrm{r}}\Delta s_i^{\mathrm{c}\mathrm{o}\mathrm{r}}\\ & =\Delta x_i^{\mathrm{a}\mathrm{f}\mathrm{f}}\Delta s_i^{\mathrm{c}\mathrm{o}\mathrm{r}}+\Delta x_i^{\mathrm{c}\mathrm{o}\mathrm{r}}\Delta s_i^{\mathrm{a}\mathrm{f}\mathrm{f}}+\Delta x_i^{\mathrm{c}\mathrm{o}\mathrm{r}}\Delta s_i^{\mathrm{c}\mathrm{o}\mathrm{r}}.\end{array}$$尽管不是尽善尽美, 但在许多情形下, 合成步$(\Delta x^{\mathrm{a}\mathrm{f}\mathrm{f}},\Delta {\lambda }^{\mathrm{a}\mathrm{f}\mathrm{f}},\Delta s^{\mathrm{a}\mathrm{f}\mathrm{f}})+(\Delta x^{\mathrm{c}\mathrm{o}\mathrm{r}},\Delta {\lambda }^{\mathrm{c}\mathrm{o}\mathrm{r}},\Delta s^{\mathrm{c}\mathrm{o}\mathrm{r}})$会比单纯的仿射变换步在下降对偶度量时表现更好.

如同在第1节中分析的理论算法, 实用算法使用了中心化步, 其中自适应地选取了中心化参数${\sigma }_k$. 仿射变换步也可以当做选取${\sigma }_k$的依据. 粗略地说,

• 若(带上步长的)仿射变换步明显地降低了对偶度量, 就没有什么必要再去做中心化, 因此可选取较小的${\sigma }_k$.
• 反之, 若下降效果不佳, 则需要较大的${\sigma }_k$将下一迭代点移向中心路径, 从而可期待从下一点出发的迭代步可以走得更远.

具体说来, 我们要计算沿着仿射变换方向的最大可行步长:$$\begin{array}{rl}{\alpha }_{\mathrm{a}\mathrm{f}\mathrm{f}}^{\mathrm{p}\mathrm{r}\mathrm{i}}& \stackrel{def}{=}\min \left(1,\underset{i:\Delta x_i^{\mathrm{a}\mathrm{f}\mathrm{f}}<0}{\min }-\frac{x_i}{\Delta x_i^{\mathrm{a}\mathrm{f}\mathrm{f}}}\right),\\ {\alpha }_{\mathrm{a}\mathrm{f}\mathrm{f}}^{\mathrm{d}\mathrm{u}\mathrm{a}\mathrm{l}}& \stackrel{def}{=}\min \left(1,\underset{i:\Delta s_i^{\mathrm{a}\mathrm{f}\mathrm{f}}<0}{\min }-\frac{s_i}{\Delta s_i^{\mathrm{a}\mathrm{f}\mathrm{f}}}\right),\end{array}$$接着定义${\mu }_{\mathrm{a}\mathrm{f}\mathrm{f}}$为用这些步长得到的对偶度量:$${\mu }_{\mathrm{a}\mathrm{f}\mathrm{f}}=(x+{\alpha }_{\mathrm{a}\mathrm{f}\mathrm{f}}^{\mathrm{p}\mathrm{r}\mathrm{i}}\Delta x^{\mathrm{a}\mathrm{f}\mathrm{f}}{)}^T(s+{\alpha }_{\mathrm{a}\mathrm{f}\mathrm{f}}^{\mathrm{d}\mathrm{u}\mathrm{a}\mathrm{l}}\Delta s^{\mathrm{a}\mathrm{f}\mathrm{f}})/n.$$中心化参数则通过以下式子得到(尚无较为坚实的理论分析, 但实际中效用较好):$$\sigma ={\left(\frac{{\mu }_{\mathrm{a}\mathrm{f}\mathrm{f}}}{\mu }\right)}^3.$$总的来说, 搜索方向的计算需要求解两个线性系统. 首先求解修正牛顿方程得到仿射变换方向, 或预估步(predictor step). 此步也用来定义计算校正步系统的右端项, 以及计算中心化参数. 接着, 搜索方向则通过求解$$\left[\begin{array}{ccc}O& A^T& I\\ A& O& O\\ S& O& X\end{array}\right]\left[\begin{array}{c}\Delta x\\ \Delta \lambda \\ \Delta s\end{array}\right]=\left[\begin{array}{c}-r_c\\ -r_b\\ -XSe-\Delta X^{\mathrm{a}\mathrm{f}\mathrm{f}}\Delta S^{\mathrm{a}\mathrm{f}\mathrm{f}}e+\sigma \mu e\end{array}\right]$$得到. 注意到上述方程组右端聚集了预估步、校正步和中心化参数, 而系数矩阵依然相同. 因此矩阵分解只需要做一次, 求解第二个系统的边际成本相对较小.

2.2 步长选取

实用算法选取步长时并不强制要求中心路径邻域${\mathcal{N}}_2,{\mathcal{N}}_{-\infty }$的作用, 而是在不打破非负性的前提下分别计算$x,s$上可以达到的最大步长. 之后再取一稍微小于这一最大值的(也不大于1的)步长. 给定迭代点$(x^k,{\lambda }^k,s^k):(x^k,s^k)>0$和迭代步$(\Delta x^k,\Delta {\lambda }^k,\Delta s^k)$, 可定义$${\alpha }_{k,\max }^{\mathrm{p}\mathrm{r}\mathrm{i}}\stackrel{def}{=}\underset{i:\Delta x_i^k<0}{\min }-\frac{x_i^k}{\Delta x_i^k},{\alpha }_{k,\max }^{\mathrm{d}\mathrm{u}\mathrm{a}\mathrm{l}}\stackrel{def}{=}\underset{i:\Delta s_i^k<0}{\min }-\frac{s_i^k}{\Delta s_i^k}.$$实用的算法于是选取位于由这些最值框定的开区间中的某个值为步长:$${\alpha }_k^{\mathrm{p}\mathrm{r}\mathrm{i}}\in (0,{\alpha }_{k,\max }^{\mathrm{p}\mathrm{r}\mathrm{i}}),{\alpha }_k^{\mathrm{d}\mathrm{u}\mathrm{a}\mathrm{l}}\in (0,{\alpha }_{k,\max }^{\mathrm{d}\mathrm{u}\mathrm{a}\mathrm{l}}),$$接着新的迭代点就是$$x^{k+1}=x^k+{\alpha }_k^{\mathrm{p}\mathrm{r}\mathrm{i}}\Delta x^k,({\lambda }^{k+1},s^{k+1})=({\lambda }^k,s^k)+{\alpha }_k^{\mathrm{d}\mathrm{u}\mathrm{a}\mathrm{l}}(\Delta {\lambda }^k,\Delta s^k).$$若$(\Delta x^k,\Delta {\lambda }^k,\Delta s^k)$修正了KKT条件前两条的不可行性, 即$$A\Delta x^k=-r_b^k=-(A{x}^k-b),A^T\Delta {\lambda }^k+\Delta s^k=-r_c^k=-(A^T{\lambda }^k+s^k-c),$$则易得在新迭代点上的不可行大小为$$r_b^{k+1}=(1-{\alpha }_k^{\mathrm{p}\mathrm{r}\mathrm{i}})r_b^k,r_c^{k+1}=(1-{\alpha }_k^{\mathrm{d}\mathrm{u}\mathrm{a}\mathrm{l}})r_c^k.$$下面的公式被许多实用算法所采纳:$${\alpha }_k^{\mathrm{p}\mathrm{r}\mathrm{i}}=\min (1,{\eta }_k{\alpha }_{k,\max }^{\mathrm{p}\mathrm{r}\mathrm{i}}),{\alpha }_k^{\mathrm{d}\mathrm{u}\mathrm{a}\mathrm{l}}=\min (1,{\eta }_k{\alpha }_{k,\max }^{\mathrm{d}\mathrm{u}\mathrm{a}\mathrm{l}}),$$这里${\eta }_k\in [0.9,1.0)$在迭代点接近原始-对偶解时将趋近于1, 以加速渐进收敛.

2.3 初始点选取

初始点的选取将极大地影响算法的强健性. 仅满足$x^0,s^0>0$随意选取$(x^0,{\lambda }^0,s^0)$往往使得算法无法收敛. 下面我们介绍一种启发式选取方法. 此法选取合理满足原始和对偶问题等式约束的初始点, 同时保证$x,s$的正性, 避免分量出现过大数值.

首先, 求得具有最小范数的满足原始约束$Ax=b$的$\tilde{x}$, 以及满足对偶约束$A^T\lambda +s=c$的$(\tilde{\lambda },\tilde{s})$, 其中$\tilde{s}$具有最小范数. 也即求解以下两个问题:$$\begin{gathered} \underset{x}{\min }\frac{1}{2}{x}^{T}x,\mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}Ax=b, \\ \underset{(\lambda ,s)}{\min }\frac{1}{2}{s}^{T}s,\mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}{A}^T\lambda +s=c. \end{gathered}$$对$A$做奇异值分解易知$\tilde{x},(\tilde{\lambda },\tilde{s})$具有显式表达:$$\tilde{x}=A^T(A{A}^T{)}^{-1}b,\tilde{\lambda }=(A{A}^T{)}^{-1}Ac,\tilde{s}=c-A^T\tilde{\lambda }.$$一般$\tilde{x},\tilde{s}$会有非正分量, 因此不适合作初始点. 我们定义$${\delta }_x=\max (-(3/2)\underset{i}{\min }{\tilde{x}}_i,0),{\delta }_s=\max (-(3/2)\underset{i}{\min }{\tilde{s}}_i,0),$$并作调整$$\hat{x}=\tilde{x}+{\delta }_{x}e,\hat{s}=\tilde{s}+{\delta }_{s}e.$$显然此时$\hat{x},\hat{s}\ge 0$. 为保证$x^0,s^0$不会太接近0且不过于相似, 我们加上两个标量$${\hat{\delta }}_x=\frac{1}{2}\frac{{\hat{x}}^T\hat{s}}{e^T\hat{s}},{\hat{\delta }}_s=\frac{1}{2}\frac{{\hat{x}}^T\hat{s}}{e^T\hat{x}}.$$注意${\hat{\delta }}_x$为$\hat{x}$的每个分量以$\hat{s}$的分量为权重时的加权平均; 类似也有${\hat{\delta }}_s$. 最后定义初始点为:$$x^0=\hat{x}+{\hat{\delta }}_{x}e,{\lambda }^0=\tilde{\lambda },s^0=\hat{s}+{\hat{\delta }}_{s}e.$$以这种方式得到初始点的计算成本大约与原始-对偶算法的一步相当.

有时, 我们有关于解的先验信息, 例如与某个相似的线性系统的解具有相同形式. 关于利用这样的"热启动(warm-start)"信息构造初始点我们放到第4节讨论.

2.4 求解线性系统

原始-对偶算法的大多计算成本在于求解线性系统. 这些系统的系数矩阵通常又大又稀疏. 这一特殊结构允许我们将它们重构成更加紧致的对称系数矩阵, 这样分解起来就更加容易便捷.

算法所要求解线性系统的一般形式如下:$$\left[\begin{array}{ccc}O& A^T& I\\ A& O& O\\ S& O& X\end{array}\right]\left[\begin{array}{c}\Delta x\\ \Delta \lambda \\ \Delta s\end{array}\right]=\left[\begin{array}{c}-r_c\\ -r_b\\ -r_{xs}\end{array}\right].$$下面对其重构. 由于$x,s$严格正, 因此对角阵$X,S$非奇异. 消去$\Delta s$并以$-X^{-1}$左乘第三行加到第一行可得$$\begin{array}{rl}\left[\begin{array}{cc}-D^{-2}& A^T\\ A& O\end{array}\right]\left[\begin{array}{c}\Delta x\\ \Delta \lambda \end{array}\right]& =\left[\begin{array}{c}-r_c+X^{-1}{r}_{xs}\\ -r_b\end{array}\right],\\ \Delta s& =-X^{-1}{r}_{xs}-X^{-1}S\Delta x,\end{array}$$其中$$D=S^{-1/2}{X}^{1/2}.$$这种形式的方程组通常被称为增广系统(augmented system). 进一步可以消去$\Delta x$, 第一行左乘$A{D}^2$加到第二行可得$$\begin{array}{rl}A{D}^2{A}^T\Delta \lambda & =-r_b-AX{S}^{-1}{r}_c+A{S}^{-1}{r}_{xs},\\ \Delta s& =-r_c-A^T\Delta \lambda ,\\ \Delta x& =-S^{-1}{r}_{xs}-X{S}^{-1}\Delta s,\end{array}$$其中$\Delta \lambda$所满足的形式被称作是正规方程组形式(normal-equations form), 这是因为它可被视作是系数矩阵为$D{A}^T$的某个线性最小二乘问题的正规方程.

大多原始-对偶算法的实施均基于上面这三个式子. 它们直接使用Cholesky分解矩阵$A{D}^2{A}^T$, 再用三角代换求解$\Delta \lambda$. $\Delta s,\Delta x$则可顺次得到. 一般的稀疏Cholesky分解软件即可应用于$A{D}^2{A}^T$, 不过要注意由于$A{D}^2{A}^T$可能会变得比较病态或奇异(这在原始-对偶算法的最后阶段尤其显著), 因此要做一些修正. Cholesky分解时可能会碰到对角线元素过小、为0甚至(由于舍入误差的存在)为负的情形. 一种解决方式是跳过这一步分解, 直接令对应$\Delta \lambda$中的元素为0.

使用正规方程的方法求解的一个缺陷是, 若$A$含稠密列, 则整个$A{D}^2{A}^T$就也是稠密的. 因此, 实用软件会识别出稠密或近似稠密的列, 将它们排除出矩阵乘积$A{D}^2{A}^T$, 然后对剩下的稀疏阵做Choleksy分解. 之后再利用诸如Sherman-Morrison-Woodbury公式的工具恢复之前排除出去的列的信息.

正规方程形式比之前两行形式的重构问题更受欢迎, 主要是因为分解稀疏对称不定矩阵的算法和软件要更加复杂, 算得更慢, 况且稀疏Cholesky分解更常用. 尽管如此, 后者仍更加简洁和灵活. 它往往可以避免由于$A$中某列稠密而导致$A{D}^2{A}^T$稠密的问题. 另外, 它允许自由变量存在. 而正规方程形式则必须增加其他的条件来表示这种变量.

2.5 实用算法

算法3 (Predictor-Corrector Algorithm)
Calculate $(x^0,{\lambda }^0,s^0)$ as described above;
for $k=0,1,2,\dots$
\quad\quad Set $(x,\lambda ,s)=(x^k,{\lambda }^k,s^k)$ and solve the primal equations for $(\Delta x^{\mathrm{a}\mathrm{f}\mathrm{f}},\Delta {\lambda }^{\mathrm{a}\mathrm{f}\mathrm{f}},\Delta s^{\mathrm{a}\mathrm{f}\mathrm{f}})$;
\quad\quad Calculate ${\alpha }_{\mathrm{a}\mathrm{f}\mathrm{f}}^{\mathrm{p}\mathrm{r}\mathrm{i}},{\alpha }_{\mathrm{a}\mathrm{f}\mathrm{f}}^{\mathrm{d}\mathrm{u}\mathrm{a}\mathrm{l}},{\mu }_{\mathrm{a}\mathrm{f}\mathrm{f}}$;
\quad\quad Set centering parameter to $\sigma =({\mu }_{\mathrm{a}\mathrm{f}\mathrm{f}}/\mu {)}^3$;
\quad\quad Solve the modified equations for $(\Delta x,\Delta \lambda ,\Delta s)$;
\quad\quad Calculate ${\alpha }_k^{\mathrm{p}\mathrm{r}\mathrm{i}},{\alpha }_k^{\mathrm{d}\mathrm{u}\mathrm{a}\mathrm{l}}$;
\quad\quad Set$$\begin{array}{rl}{x}^{k+1}& =x^k+{\alpha }_k^{\mathrm{p}\mathrm{r}\mathrm{i}}\Delta x,\\ ({\lambda }^{k+1},s^{k+1}& =({\lambda }^k,s^k)+{\alpha }_k^{\mathrm{d}\mathrm{u}\mathrm{a}\mathrm{l}}(\Delta \lambda ,\Delta s);\end{array}$$end (for)

对于上述算法, 尚无可用的收敛理论. 事实上, 存在使上述预估-校正算法失效的例子. 尽管我们可以加上一些防护措施将它强制嵌入现已存在的收敛算法或者改善它的强健性, 但因失效比较少见, 许多实用代码并不关心这一点.

当给定的线性规划不可行或无界, 则算法往往会发散, 其中不可行度量$r_b^k,r_c^k$或对偶度量${\mu }_k$可能会趋于$\infty$. 由于不可行和无界很好识别, 因此内点法的代码都内含启发程序探测、报告这些现象. 探测不可行和无界的更严格化的方法则会使用齐次自对偶(homogeneous self-dual formualtion).

3. 其他原始-对偶算法和相关扩展

3.1 其他路径跟踪算法

框架1是许多其他路径跟踪类算法的基础. 从实用观点上看它们可能不那么重要, 不过因其优美和强大的理论性质, 我们在此简略涉及.

1. 一些路径跟踪算法保守地选取中心化参数$\sigma$(即$\sigma$没比1小多少), 这样就使得算法在得到的方向上总是能取到接近单位的步长. 这种算法被称作是短步路径跟踪算法(short-step path-following). 正是由于它们要求迭代点一直在${\mathcal{N}}_2$里徘徊, 它们趋向解的进程也就较为缓慢. 不过从理论的角度, 它们具有更好的复杂度估计(对此有类似定理4的结论, 不过$n$要用$n^{1/2}$替换).
2. 对使用两个${\mathcal{N}}_2$邻域的预估-校正算法我们有更好的结论. 这两个邻域中有一个是嵌在另一个里头的. 尽管都叫预估-校正, 但此"预估-校正"非彼预估-校正. 此算法每隔一步就是预估步, 从内部的邻域中某点出发, 沿着仿射变换方向移动(即令$\sigma =0$)至外邻域的边界. 邻域边界之间的间隔使这一步能带来$\mu$上充分的下降. 与预估步交替进行的是校正步(令$\sigma =1,\alpha =1$), 它将下一个迭代点拉回到内部邻域中, 为下一次预估步做准备. 此预估-校正算法产生的对偶度量序列${\mu }_k$超线性收敛到0, 而大多其他算法仅有线性收敛.

3.2 势能下降算法

势能下降算法potential-reduction methods所取的迭代步形式类似于路径跟踪算法, 但它们并不显式地跟踪中心路径$\mathcal{C}$. 事实上它们与中心路径是独立的. 它们使用对数势能函数来度量${\mathcal{F}}^o$中每个点的价值并期望每步迭代后均能得到在这个函数上的一定量的下降. 我们以$\Phi$表示原始-对偶势能函数, 它通常具有两个重要性质:$$\begin{array}{rl}\Phi & \to \infty \Leftarrow \exists i:x_i{s}_i\to 0;\mu =x^{T}s/n̸=0,\\ \Phi & \to -\infty \iff (x,\lambda ,s)\to \Omega .\end{array}$$第一条性质防止任何一对乘积$x_i{s}_i$独立地趋向于0, 从而不让迭代点聚于非负象限边界上. 第二条性质将$\Phi$与解集$\Omega$联系起来. 若算法强制$\Phi \to -\infty$, 则此性质就保证迭代序列是趋近于解集中的点的.

一种有趣的原始-对偶势能函数定义为$${\Phi }_{\rho }(x,s)=\rho \log x^{T}s-\sum _{i=1}^n\log x_i{s}_i,$$其中参数$\rho >n$. 诸如所有基于框架1的算法, 势能下降算法对某个${\sigma }_k\in (0,1)$求解修正牛顿方程, 再沿着得到的方向选取步长${\alpha }_k$. 例如, 我们可以选取步长${\alpha }_k$为沿着搜索方向近似极小化${\Phi }_{\rho }$的值. 若固定${\sigma }_k\equiv n/(n+\sqrt{n})$, 我们可保证每一步在${\Phi }_{\rho }$上均有固定量的下降. 这样${\Phi }_{\rho }\to -\infty$, 得收敛. 匹配保守选取方式的自适应和启发式的${\sigma }_k,{\alpha }_k$选取也在我们的考虑范围内.

3.3 相关扩展

线性规划的原始-对偶算法可以扩展至应用于更多的问题. 算法可以简单扩展至求解单调线性互补问题(monotone linear complementarity problem, monotone LCP)和凸二次规划问题, 其中保留了线性规划时的收敛性和多项式复杂度.

3.3.1 单调线性互补问题

单调LCP是要求${\mathbb{R}}^n$中的$x,s$满足以下条件:$$s=Mx+q,(x,s)\ge 0,x^{T}s=0,$$其中$M$为$n\times n$的半正定矩阵, $q\in {\mathbb{R}}^n$. 上面条件与KKT条件显然很相似: 后两个条件分别对应于KKT中的非负性条件与互补松弛条件, 而$s=Mx+q$则类似于稳定性条件与可行性条件.

求解单调LCP的内点法与线性规划算法联系紧密. 对偶度量被重新定义为互补度量(complementarity measure), 定义式仍然是$\mu =x^{T}s/n$. 上述条件可以被类似地用对角阵写作$$\left[\begin{array}{c}Mx+q-s\\ XSe\end{array}\right]=0,(x,s)\ge 0.$$对此路径跟踪步的一般形式亦可类似地定义为:$$\left[\begin{array}{cc}M& -I\\ S& X\end{array}\right]\left[\begin{array}{c}\Delta x\\ \Delta s\end{array}\right]=\left[\begin{array}{c}-(Mx+q-s)\\ -XSe+\sigma \mu e\end{array}\right],$$其中$\sigma \in [0,1]$. 类似亦可给出形如第2节使用算法的扩展.

3.3.2 凸二次规划问题

对于凸二次规划的扩展我们放到后面讲述二次规划的章节中讨论. 而对一般非线性规划问题的扩展则是第十九章的主题.
内点法在求解半定规划问题时是极其高效的. 半定规划问题是一类涉及对称半定矩阵变量的问题. 半定规划自上世纪九十年代早期开始便成为了研究焦点, 它在许多领域包括控制论、组合优化上都有应用.

4. 展望和软件

上世纪八十年代内点法的出现挑战了单纯形法在求解线性规划问题上的统治地位. 大约九十年代前, 包含本章第二节中的内点法代码问世, 且在那时在许多大型问题上的表现要优于单纯形法. 之后单纯形法也不甘示弱, 相继出现了许多改进和软件包, 如CPLEX, XPRESS-MP. 这些改进要归功于算法上如最速下降边旋转和改良的定价方式的提出, 也离不开对算法细节上的反复推敲. 与此同时, 内点法的效率也不断提高, 这主要通过改进求解方程的算法以及使用更高阶的校正子. 这一时期涌现了许多好的内点法代码, 例如PCx, HOPDM, BPMPD, LIPSOL. 它们已被投入许多应用当中.

一般, 单纯形法在中小型问题上较快, 而内点法在较大型的问题上更具优势. 但这并不是一定的. 算法的效果相当依赖于问题的特定结构. 内点法往往不能利用好解的先验信息, 例如对解本身的估计或最优基的估计. 因此, 在"热启动"存在时, 内点法可能比单纯形法表现得要差. 此类情形的一个例子是求解整数规划的分支定界法, 其中确定分支定界树的每个节点需要父节点的信息. 此外, 我们可能希望求解一系列线性规划问题, 其中我们对数据作轻微扰动以观察解对扰动的灵敏度, 或者是我们通过一系列线性规划来逼近一个非线性规划问题. Yildirim和Wright提出如何改进给定点(如问题的近似解)以获取理论上合理的初始点, 其中伴随根据给定点质量证明的复杂度理论. 但实际上, 这些技术仅能带来算法表现上的微小改进.

内点法的好处之一是它(相比于单纯形法)易于编程实现. 这里面最复杂的操作是求解大型线性系统, 然而这在许多软件中是可直接调用的.