你怎么不牵我！

# 线性规划问题 (Linear Programming, LP )

线性规划问题的规范形式：

$$\begin{cases} min & \textbf{c}^T\textbf{x}\\ s.t. & \textbf{Ax}\geq\textbf{b}\\ & \textbf{x}\geq\textbf{0} \end{cases}$$

线性规划问题的标准形式：

$$\begin{cases} min & \textbf{c}^T\textbf{x}\\ s.t. & \textbf{Ax}=\textbf{b}\\ & \textbf{x}\geq\textbf{0} \end{cases}$$

上述两种形式可以相互转化。譬如$x_1+x_2\geq 3$ 可以通过引入一个新的非负变量$x_3\geq 0$ 使得变为$x_1+x_2-x_3= 3$。

我们考虑标准形式的 LP，假设$rank(A)=m$（$A$ 是$m\times n$ 的），故$A$ 中有$m$ 个线性无关的列向量，它们构成了满秩方阵$B_{m\times m}$。再把$A$ 中其余各列
组成的子阵记为$N$，即$A=(B,N)$，再把$\textbf{x}=(x_1,x_2,...,x_n)^T$ 分为两部分：记为$\textbf{x}_B$ 和$\textbf{x}_N$，则有：

$$B\textbf{x}_B+N\textbf{x}_N=\textbf{b}$$

可以解得

$$\textbf{x}_B=B^{-1}\textbf{b}-B^{-1}N\textbf{x}_N$$

那么我们可以把$\textbf{x}_N$ 当作是自由变量，可以取任意值$\bar{\textbf{x}}_N$，然后代入得到$\bar{\textbf{x}}_B$ 后就得到了原方程的一组解。

证明：若$\textbf{x}$ 是基本可行解，那么它只有非负分量而且对应的向量为基向量，故线性无关。

反之，由于$\textbf{x}$ 是可行解，故有$\textbf{x}\geq 0$ 且$A\textbf{x}=\textbf{b}$。不失一般性假设$\textbf{x}$ 的前$k$ 个分量为正：

$$\textbf{x}=(\bar{x}_1,\bar{x}_2,...,\bar{x}_k,0,...,0)^T\\ \bar{x_i}>0,i=1,...,k$$

那么显然有$k\leq m$，可以讨论：

• 若$k=m$，则自然$\textbf{x}$ 就是对应于$B=(a_1,a_2,...,a_k)$ 的一个基本可行解，其中$a_i$ 就是分量$\bar{x}_i$ 对应的$A$ 中的列向量。
• 若$k<m$，此时需要人为再添加几个列向量。因为$rank(A)=m$，那么$A$ 中至少有$m$ 个线性无关的列向量，故可以从$A$ 中选取$m-k$ 个线性无关的列向量，与之前的$k$ 个列向量组成了$B$。此时，$\textbf{x}_B=(\bar{x}_1,...,\bar{x}_k,0,...,0)^T$ 是基本可行解。

---

对于一个基本可行解$\bar{\textbf{x}}$，如果它的所有基分量都取正值，那么称它为非退化的，否则如果有的基分量是 0，就称它为退化的。
一个可行基对应一个基本可行解，反之若一个基本可行解是非退化的，那么也唯一对应着一个可行基。

一般地说，如果一个基本可行解是退化的，那么它可由不止一个可行基得到。一个标准形式的 LP 问题最多有$\begin{pmatrix}n\\m\end{pmatrix}$ 个可行基，从而基本可行解的个数不会超过$\begin{pmatrix}n\\m\end{pmatrix}$ 个，解集空间的凸多面体顶点数也不超过$\begin{pmatrix}n\\m\end{pmatrix}$ 个。一个 LP 问题，如果它的所有基本可行解都是非退化的，就说该问题是非退化的，否则称之为退化的。

证：设$\textbf{x}^0=(x_1^0,...,x_n^0)^T$ 是一个可行解，则有$A\textbf{x}_0=\textbf{b},\textbf{x}^0\geq 0$。不失一般性，设$\textbf{x}^0$ 的正分量为前$k$ 个。如果它们对应的$A$ 中的列向量$A_1,...,A_k$ 已经线性无关了，则根据前述定理，$\textbf{x}^0$ 已经是一个基本可行解了。

否则，若不线性无关，那么存在不全为零的数$\delta_j,j=1,...,k$ 使得：

$$\sum_{j=1}^k\delta_jA_j=\textbf{0}$$

令$\delta_j=0,j=k+1,...,n$ 得到向量$\delta=(\delta_1,...,\delta_k,0,...,0)^T$, 故有：

$$A\delta=\textbf{0}$$

由于$x_j^0>0,j=1,...,k$，故可以选择一个足够小的正数$\varepsilon$ 使得：

$$\textbf{x}^0\pm\varepsilon\delta\geq 0$$

故有

$$A(\textbf{x}^0\pm\varepsilon\delta)=A\textbf{x}^0+\varepsilon A\delta=\textbf{b}$$

而我们可以通过选择$\varepsilon$，使得$\textbf{x}^0+\varepsilon\delta$ 和$\textbf{x}^0-\varepsilon\delta$ 中，有一个可以比$\textbf{x}^0$ 的非零分量少一个。故一旦正分量对应的$A$ 中列向量不线性无关，我们就可以重复这个过程，直到非零分量剩下一个时就肯定是线性无关的了（或在某一次操作后就线性无关了）。证毕。

证：设$\textbf{x}^0$ 是一个最优解，则可以根据上个证明中的过程构造出两个可行解，并且它们的目标函数值为：

$$\textbf{c}^T(\textbf{x}^0-\varepsilon\delta)=\textbf{c}^T\textbf{x}^0-\varepsilon\textbf{c}^T\delta\\ \textbf{c}^T(\textbf{x}^0+\varepsilon\delta)=\textbf{c}^T\textbf{x}^0+\varepsilon\textbf{c}^T\delta$$

因为$\textbf{c}^T\textbf{x}^0$ 是最优解，故有：

$$\textbf{c}^T(\textbf{x}^0\pm\varepsilon\delta)\geq\textbf{c}^T\textbf{x}^0\\$$

即$\varepsilon\textbf{c}^T\delta=0$，故得到的两个新的可行解$\textbf{x}^0\pm\varepsilon\delta$ 的目标函数值也是最优的。同理可以取$\varepsilon$ 使得其中非零分量减少一个，重复操作直到其成为一个基本可行解。证毕。

# LP 问题的单纯形法

既然在上一部分已经证明了最优解都可以转化为基本可行解，故我们只需要讨论基本可行解就足够了，不会漏掉最优情况。
这里有两个问题需要解决：一是如何找到一个基本可行解，二是如何判断基本可行解是否最优了，以及如何迭代转化基本可行解。前一个问题留在下一章解决。

假设已经找到了一个非退化的基本可行解$\bar{\textbf{x}}$，即找到了一个基$B$，此时可以将方程组$A\textbf{x}=\textbf{b}$ 转化为与之等价的方程组：

$$\textbf{x}_B+B^{-1}N\textbf{x}_N=B^{-1}\textbf{b}$$

不失一般性假定$B=(A_1,...,A_m)$，记向量：

$$\bar{A_j}=B^{-1}A_j\\ \bar{\textbf{b}}=B^{-1}\textbf{b}$$

则可以把方程组转化为：

$$\textbf{x}_B+B^{-1}N\textbf{x}_N=\bar{\textbf{b}}$$

这个方程组被称为对应于基$B$ 的典则方程组，简称典式。
相应的，可以改写目标函数值为：

$$\begin{aligned} \textbf{c}^T\textbf{x} &=\textbf{c}_B^T\textbf{x}_B+\textbf{c}_N^T\textbf{x}_N\\ &=\textbf{c}_B^T(\bar{\textbf{b}}-B^{-1}N\textbf{x}_N)+\textbf{c}_N^T\textbf{x}_N\\ &=\textbf{c}_B^T\bar{\textbf{b}}-\sum_{j=m+1}^n(\textbf{c}_B^T\bar{A}_j-c_j)x_j \end{aligned}$$

其中不失一般性地假设了$\textbf{x}$ 的前$m$ 个分量对应$\textbf{x}_B$，后$n-m$ 个分量对应了$\textbf{x}_N$。可以注意到：

$$(\bar{A}_1,...,\bar{A}_m)=B^{-1}(A_1,...,A_m)=B^{-1}B=I$$

故$\bar{A}_i$ 其实就是一个只有第$i$ 个分量为 1，其他分量都为 0 的列向量。故有：

$$\textbf{c}_B^T\bar{A}_j-c_j=0,j=1,...,m$$

引入记号：

$$\zeta_j=\textbf{c}_B^T\bar{A}_j-c_j,j=1,...,n$$

向量表示为:

$$\zeta^T=\textbf{c}_B^TB^{-1}A-\textbf{c}^T=(\zeta_B,\zeta_N)^T\\ =(0,\zeta_N)^T$$

从而目标函数值为：

$$\textbf{c}^T\textbf{x}=\textbf{c}_B^T\bar{\textbf{b}}-\zeta^T\textbf{x}$$

证：对于任意一个可行解$\textbf{x}$，一定有$\textbf{x}\geq 0$，那么就有：

$$\textbf{c}^T\bar{\textbf{x}}=\textbf{c}_B^T\bar{\textbf{b}}\leq\textbf{c}_B^T\bar{\textbf{b}}-\zeta^T\textbf{x}=\textbf{c}^T\textbf{x}$$

证：令$\textbf{d}=\begin{pmatrix}-\bar{A}_k\\0\end{pmatrix}+e_k$，其中$e_k$ 是第$k$ 个分量为 1，其余分量为 0 的列向量。
因为$\bar{A}_k\leq 0$，所以有$\textbf{d}\geq 0$。而：

$$\begin{aligned} A\textbf{d} &=(B,N)\begin{pmatrix}-\bar{A}_k\\0\end{pmatrix}+Ae_k\\ &=-A_k+A_k\\ &=\textbf{0} \end{aligned}$$

对于充分大的数$\theta$，观察$\bar{\textbf{x}}+\theta\textbf{d}$，此时有：

$$A(\bar{\textbf{x}}+\theta\textbf{d})=A\bar{\textbf{x}}+\theta A\textbf{d}=\textbf{b}\\ \bar{\textbf{x}}+\theta\textbf{d}\geq 0$$

所以$\bar{\textbf{x}}+\theta\textbf{d}$ 是原问题的可行解，那么目标函数值为：

$$\begin{aligned} \textbf{c}^T(\bar{\textbf{x}}+\theta\textbf{d})&=\textbf{c}^T\bar{\textbf{x}}+\theta(\textbf{c}_B^T,\textbf{c}_N^T)\begin{pmatrix} -\bar{A}_k\\ 0 \end{pmatrix}+\theta\textbf{c}^Te_k\\ &=\textbf{c}^T\bar{\textbf{x}}-\theta(\textbf{c}_B^T\bar{A}_k-c_k)\\ &=\textbf{c}^T\bar{\textbf{x}}-\theta\zeta_k \end{aligned}$$

又因为$\zeta_k>0$，故可以无限大地选择$\theta$，使得目标函数无限小，故无界，证毕。

证：只需要找出$\hat{\textbf{x}}$ 是什么。
类似前述证明，令

$$\textbf{d}=\begin{pmatrix}-\bar{A}_k\\0\end{pmatrix}+e_k$$

则有$A\textbf{d}=\textbf{0}$。令：

$$\hat{\textbf{x}}=\bar{\textbf{x}}+\theta\textbf{d}=\begin{pmatrix}\bar{\textbf{b}}-\theta\bar{A}_k\\0\end{pmatrix}+\theta e_k$$

我们选取合适的$\theta$。首先为了$\hat{\textbf{x}}\geq 0$, 可以选

$$\theta=min\{\frac{\bar{b}_i}{\bar{a}_{ik}}|\bar{a}_{ik}>0,i=1,...,m\}\\ =\frac{\bar{b}_r}{\bar{a}_{rk}}$$

从而类似前述证明，$\hat{\textbf{x}}$ 是一个可行解。
先考虑$\hat{\textbf{x}}$ 的各个分量：

$$\begin{aligned} & \hat{x}_i=\bar{b}_i-\frac{\bar{b}_r}{\bar{a}_{rk}}\bar{a}_{ik},i=1,...,m,i\neq r\\ & \hat{x}_r=0\\ & \hat{x}_k=\theta=\frac{\bar{b}_r}{\bar{a}_{rk}}\\ & \hat{x}_j=0,j=1,...,n,j\neq k \end{aligned}$$

为证明其是一个基本解，则需要证明

$$A_1,...,A_{r-1},A_k,A_{r+1},...,A_m$$

是线性无关的。（即把第$r$ 个分量换出去，把第$k$ 个换进来）。
反设不是线性无关的，由于先前的$A_1,...,A_r,...,A_m$ 是线性无关的，那么就会有：

$$A_k=\sum_{i=1,i\neq r}^my_iA_i$$

同时有：

$$A_k=B\bar{A}_k=\sum_{i=1}^m\bar{a}_{ik}A_i$$

将上述两式相减可得：

$$\bar{a}_{rk}A_r+\sum_{i=1,i\neq r}^m(\bar{a}_{ik}-y_i)A_i=0$$

已知有$\bar{a}_{rk}> 0$，故与 “$A_1,...,A_r,...,A_m$ 线性独立” 矛盾了。
最后，由非退化假设$\bar{b}>0$，因而$\theta>0$，故：

$$\textbf{c}^T\hat{\textbf{x}}=\textbf{c}^T\bar{\textbf{x}}-\theta\zeta_k<\textbf{c}^T\bar{\textbf{x}}$$

证毕。

由以上定理我们知道相应于基本可行解$\bar{\textbf{x}}$ 的向量$\zeta^T=\textbf{c}_B^TB^{-1}A-\textbf{c}^T$ 有重要的作用，我们称他为基本可行解$\bar{\textbf{x}}$ 的检验数向量，它的各个分量称为检验数。

上述定理给了一个从一个基本可行解$\bar{\textbf{x}}$ 变成一个更好的的基本可行解$\hat{\textbf{x}}$ 的办法，即用$x_k$ 代替原来的基变量$x_r$ 成为新的基。

当线性规划原问题是退化问题时，由线性规划问题的几何解释可知，通过该可行域某个极点的超平面超过 n 个，所以该点为一个退化的极点。根据摄动法原理，可在退化问题约束方程的右边项做微小的扰动，使得超平面有一个微小的位移，原来相交于一点的若干个超平面略微错开一些，退化极点变成不退化极点。决策者可根据问题的实际情况，适当增加或减少某些资源的数量，使得其迭代变为非退化的，以得到问题的最优解。

单纯形法步骤：

1. 找一个初始的可行基$B$。
2. 求出对应的典式及检验数向量$\zeta$。
3. 求$\zeta_k=max\{\zeta_j|j=1,...,n\}$。
4. 若$\zeta_k\leq 0$，则停止迭代，得到最优解$\textbf{x}=\begin{pmatrix}\textbf{x}_B\\ \textbf{x}_N\end{pmatrix}=\begin{pmatrix}\bar{\textbf{b}}\\0\end{pmatrix}$。
5. 若$\bar{A}_k\leq 0$ 则停止，原问题无界。
6. 求$\frac{\bar{b}_r}{\bar{a}_{rk}}=min\{\frac{\bar{b}_i}{\bar{a}_{ik}}\;|\;\bar{a}_{ik}>0,i=1,...,m\}$。
7. 用$A_k$ 代替$A_{B_r}$，得到新的基转第二步。

# LP 问题的初始解

不失一般性，可以在方程左右乘以 (-1) 使得$\textbf{b}\geq 0$，设原问题为：

$$\begin{cases} min\; &z=\textbf{c}^T\textbf{x}\\ s.t. & A\textbf{x}=\textbf{b}\quad(\textbf{b}\geq 0)\\ & \textbf{x} \geq 0 \end{cases}$$

我们要找一个初始的基本可行解，一个很自然的想法就是为啥不直接对$(A,\textbf{b})$ 行变换，然后行变换出来一个单位矩阵后，就得到了一个基本可行解。但实际上无法确认究竟行变换后$\textbf{b}$ 是否能保证不为负数。譬如我不能先入为主就认为$x_1,x_2,x_3\neq 0,x_4,x_5=0$ 是一组可行解，很可能做完行变换后$\textbf{b}$ 就不是非负的了。

为了解决 "得到一个初始的基本可行解" 这个问题，我们引入$m$ 个人工变量：

$$\textbf{x}_\alpha=(x_{n+1},...,x_{n+m})^T$$

并研究以下 LP 问题：

$$\begin{cases} min\; & g=\sum_{i=n+1}^{n+m}x_i\\ s.t. & A\textbf{x}+\textbf{x}_\alpha=\textbf{b}\\ & \textbf{x}\geq 0,\quad\textbf{x}_\alpha\geq 0 \end{cases}$$

不难发现，原问题有可行解当且仅当上述 LP 问题有最优解$g=0$。此时即$\textbf{x}_\alpha=0$。
而上述 LP 问题已经有了一个天然的初始基本可行解，即后$m$ 个变量。

但是有一种情况就是有最优解$g=0$，人工添加的$\textbf{x}_\alpha=0$，但是某个人工添加的变量$x_{n+i}$ 取值为 0，但却是基变量。

此时会发现有一行（即$x_{n+i}$ 对应的那一行）有$\bar{b}_i=0$，因为那一行的基变量的系数只有$x_{n+i}$ 是 1，其他基变量的系数都是 0。
而不是基变量的变量取值是 0，又因为最优解下$x_{n+i}$ 取值也是 0，故那一行结果$\bar{b}_i=0$。
这样的话，我们就可以考虑那一行中，如果$x_1,...,x_n$ 中有系数不是 0 的，就把他换进基，把$x_{n+i}$ 换出去。因为$\bar{b}_i=0$，可以证明这样的换基是不改变目标函数的。
如果全是 0，说明原线性约束亏秩，可以直接把这一行删除掉，是无用约束。

# LP 问题的对偶性

对偶问题其实有很直观的理解。考虑如下一个例子：

$$\begin{cases} min\;& 2x_1+3x_2\\ s.t. & x_1+x_2\geq 1\\ & x_2\geq 2\\ & x_1,x_2\geq 0 \end{cases}$$

我们很显然把第一个不等式乘以 2 加上第二个不等式就得到了$2x_1+3x_2\geq 4$。
我们就得到了一个目标函数的一个下界，而这个下界是可能取不到的，譬如这个例子就是。
但是目标函数一定是大于等于这个下界的。
所以对偶问题实际上就是在寻找一个原问题的下界。

理解一下，$A^T\textbf{y}$ 就是把各个不等式加权加起来，而且要求结果是小于等于$\textbf{c}$ 的就表示加起来后，右边大于等于的东西要严格小于等于目标函数值。然后$\textbf{b}^T\textbf{y}$ 就是加权后右边大于等于的东西，他一定是目标函数的下界。而取 max 实际上就是取下界中最大的一个。

证：考虑标准形式的 LP 问题和其对偶问题，对于任意的可行解$\textbf{x},\textbf{y}$：

$$\begin{aligned} & \because A\textbf{x}= b\\ & \therefore \textbf{y}^TA\textbf{x}=\textbf{y}^T\textbf{b}\\ & \because A^T\textbf{y}\leq \textbf{c},\;\textbf{x}\geq 0\\ & \therefore \textbf{y}^TA\textbf{x}\leq\textbf{c}^T\textbf{x}\\ & \therefore \textbf{y}^T\textbf{b}=\textbf{y}^T A\textbf{x}\leq\textbf{c}^T\textbf{x}\\ \end{aligned}$$

证明了对偶问题的目标函数总是小于原目标函数的。
现假设原问题有了一个最优解为$\hat{\textbf{x}}$，它对应的基为$\hat{B}$（即$\hat{\textbf{x}}_B=\hat{B}^{-1}\textbf{b}$），那么$\textbf{y}^T=\textbf{c}_{\hat{B}}^T\hat{B}^{-1}$ 就是其对偶问题的一个可行解：

$$A^T\textbf{y}-\textbf{c}=A^T(\textbf{c}_{\hat{B}}^T\hat{B}^{-1})^T-\textbf{c}=\zeta\leq 0\\$$

此时它们的目标函数值相同：

$$\textbf{y}^T\textbf{b}= \textbf{c}_{\hat{B}}^T\hat{B}^{-1}\textbf{b}=\textbf{c}_{\hat{B}}^T\hat{\textbf{x}}_B$$

故此时对偶问题一定是取到了可能的最上界，即为最优值。证毕。

其实原始的 LP 问题和其对偶问题只有以下几种解的情况：

1. 两个问题都有最优解，且最优解相等。
2. 一个问题无界，另一个问题无可行解。
3. 两个问题都无可行解。

证：考虑规范形式，就会有：

$$\begin{aligned} & \because\textbf{x},\textbf{y},(A\textbf{x}-\textbf{b}),(\textbf{c}-A^T\textbf{y})\geq 0\\ & \therefore \textbf{y}^T(A\textbf{x}-\textbf{b})\geq 0,\textbf{x}^T(\textbf{c}-A^T\textbf{y})\geq 0\\ & \because\textbf{y}^T(A\textbf{x}-\textbf{b})+\textbf{x}^T(\textbf{c}-A^T\textbf{y})=\textbf{x}^T\textbf{c}-\textbf{y}^T\textbf{b} \end{aligned}$$

所以$\textbf{x}^T\textbf{c}=\textbf{y}^T\textbf{b}\Leftrightarrow\textbf{y}^T(A\textbf{x}-\textbf{b})=0,\textbf{x}^T(\textbf{c}-A^T\textbf{y})=0$，证毕。

# LP 问题的计算复杂度

• 苏联数学家 Khachian 于 1979 年提出椭球算法理论上证明了 LP 是一个多项式可解问题。
• 1984 年，N.Karmarkar 提出了 Karmarkar算法 。单纯形法是从可行域的一个顶点开始，沿着可行域的边从一个顶点转移到另一个顶点。但是 Karmarkar算法 却是从多面体内部一个点开始，产生一个直接穿过多面体内部的点列：
  
  Karmarkar 证明了该算法迭代次数的上届是$O(nL)$ 每次迭代的平均计算量为$O(n^{2.5})$，因而其算法的总的时间复杂度为$O(n^{3.5}L)$，其中$n$ 是维数，$L$ 为输入规模。
• 之后也诞生了诸如仿射尺度法、路径追踪法等 LP 问题的解法，其中最好的时间复杂度在$O(n^3L)$。已经在商业化软件中实现，但是内点法有一个固有弱点就是仅提供近似最优解。需另外加一个纯化的过程。
• 单纯形法有 “热启动” 的特性，即从最后求得的基再开始能大大缩短求解时间。

# 整数线性规划

# 整数线性规划问题 ILP

当我们想求解一个 ILP 问题时，可以想象先求解对应的 LP 问题，然后自然地把结果舍入到一个整数解。但是从 LP 问题的可行解舍入到一个可行的整数解是困难的，甚至是不可行的：

实际上这个舍入的过程复杂度已经相当于求解原问题。

# Gomory割平面法

考虑一种特殊的整数线性规划问题，即纯整数线性规划问题 ($P$)：

$$\begin{cases} min\;& z=\mathbf{c}^T\mathbf{x}\\ s.t.\;& A\textbf{x}=\textbf{b}\\ & \textbf{x}\geq 0\\ \end{cases}$$

其中，$A,\mathbf{c},\mathbf{x},\mathbf{b}$ 全部为整数向量。（其实这里并不失一般性，因为显然可以同乘$A,\mathbf{b}$ 的分母最小公倍数使其为整数，而$\mathbf{c}$ 怎么乘也不会影响最优解。）

对于这类问题，我们不妨记问题$P$ 的可行解范围为$D$, 当$D\neq\emptyset$ 时，它是由有限个或可数个格点构成的集合。

现在我们去掉 “$\mathbf{x}$ 是整数向量这个条件”，得到了一个 LP 问题$(P_0)$ 称为$P$ 的松弛问题，记$P_0$ 的可行域为$D_0$，很显然有以下性质：

1. $D\subset D_0$。
2. 若$P_0$ 无可行解，则$P$ 无可行解。
3. $P_0$ 的最优值是$P$ 最优值的一个下界。
4. 若$P_0$ 的最优解是一个整数向量，则他也是$P$ 的最优解。

Gomory割平面法 的基本思想就是，先用单纯形法求解$P_0$，然后若最优解$\mathbf{x}^*$
不是整数向量，那么就添加一个约束条件（称为割平面条件），在保证不割去格点（$D$ 保持不变）的前提下，割去$\mathbf{x}^*$，然后再求解$P_0'$，如此循环直到最优解是整数向量为止。

如果在增加约束的过程中，得到的 LP 问题无可行解，则原 ILP 问题也无可行解。如果 LP 问题无解，则预案 ILP 问题无可行解或无界。

---

下面我们给出生产这样的割平面条件的代数方法。

给定一个$P_0$ 的最优基本可行解$\mathbf{x}^0$，设它对应的基为$B=(A_{B_1},...,A_{B_m})$，而$x_{B_1},...,x_{B_m}$ 是最优基本可行解的基变量，基变量的下标集合为$S$，非基变量的下标集合为$\bar{S}$。那么最优解的典式满足：

$$x_{B_i}+\sum_{j\in\bar{S}}\bar{a}_{ij}x_j=\bar{b_i},i=1,...,m$$

若$\bar{b}_i,i=1,...,m$ 全为整数时，则已经找到了一个 ILP 的最优解。否则至少有一个$\bar{b}_l$ 不为整数。

由于上述典式中，$x_{B_i},x_j$ 都是非负的，因此有：

$$\sum_{j\in\bar{S}}[\bar{a}_{ij}]x_j\leq\sum_{j\in\bar{S}}\bar{a}_{ij}x_j\\ x_{B_i}+\sum_{j\in\bar{S}}[\bar{a}_{ij}]x_j\leq\bar{b_i},i=1,...,m$$

然后我们想保证原 ILP 问题的可行域不变。很显然对于上述下取整后的式子，由于左侧都是整数，那么我们对右侧也下取整，并不影响不等式：

$$x_{B_i}+\sum_{j\in\bar{S}}[\bar{a}_{ij}]x_j\leq[\bar{b}_i],i=1,...,m$$

将原来的等式和上述不等式相减，得到：

$$\sum_{j\in\bar{S}}(\bar{a}_{lj}-[\bar{a}_{lj}])x_j\geq \bar{b}_l-[\bar{b}_l]$$

于是我们就得到了一个新的形如$fx\geq b$ 的线性约束，称它为对应于生成行$l$ 的 Gomory 割平面条件。
由于它是一个大于等于的式子，我们还需要引入新的非负松弛变量$s$，使得：

$$-\sum_{j\in\bar{S}}f_{lj}x_j+s=-f_l$$

证明：原问题可行域不变在于对$\bar{b}_l$ 下取整的那一步，原问题显然不会受到影响。但对于新的松弛问题，显然有$-f_l<0$ (因为$\bar{b}_l$ 不是整数)，那么由于引入了松弛变量$s$，新的基本解就成了：

$$\begin{aligned} \mathbf{x}'^T&=(x_{B_1},...,x_{B_m},s,x_{\bar{S}_1},...,x_{\bar{S}_n})\\ &=(\bar{b}_1,...,\bar{b}_m,-f_l,0,...,0) \end{aligned}$$

显然不再可行，因为中间有一个$-f_l$ 负数。

---

因此我们可以总结割平面法的求解过程：

1. 舍去$\mathbf{x}$ 为整数向量的要求，给出松弛问题。
2. 使用单纯形法求出松弛问题的一个最优解$\mathbf{x}^0$。
3. 若第$\mathbf{x}_l^0$ 不为整数，则引入新的非负松弛变量$s$，并添加约束$-\sum_{j\in\bar{S}}f_{lj}x_j+s=-f_l$。
4. 继续求解松弛问题，直到$\mathbf{x}^0$ 为整数向量。

# 分枝定界法

还是考虑纯整数规划问题 ($P$)：

$$\begin{cases} min\;& z=\mathbf{c}^T\mathbf{x}\\ s.t.\;& A\textbf{x}=\textbf{b}\\ & \textbf{x}\geq 0\\ \end{cases}$$

其中，$A,\mathbf{c},\mathbf{x},\mathbf{b}$ 全部为整数向量。同样是求解松弛问题$P_0$，若最优解中有一维不为整数，即$x_l^*$ 不为整数的话，我们可以考虑将松弛问题划分为两个子问题：

$$\begin{cases} min\;& z=\mathbf{c}^T\mathbf{x}\\ s.t.\;& A\textbf{x}=\textbf{b}\\ & \textbf{x}\geq 0\\ & x_l\leq [x_l^*]\\ \end{cases} \quad \begin{cases} min\;& z=\mathbf{c}^T\mathbf{x}\\ s.t.\;& A\textbf{x}=\textbf{b}\\ & \textbf{x}\geq 0\\ & x_l\geq [x_l^*]+1\\ \end{cases}$$

显然这两个问题是互斥的，可行域不相交。因此我们可以分别求解这两个子问题，然后取其中最优解即可。于是会有分枝树：

如果原问题有界，则显然分枝的过程不会无限进行下去。分枝会因为两种情况而停下来：得到了一组整数最优解，或松弛问题不可行。这样最后把所有叶子结点取个最优值即得到原问题的最优值。

很显然我们也可以进行剪枝优化，如果目前分枝树中已经得到了整数最优解，且最优值为$z^*$，那么对于一个非叶子结点，如果它的松弛问题的最优值已经不优于$z^*$ 了就可以停止对它分枝，称它为死点。

# 非线性规划

# 基本概念

# 凸函数和凸规划

很显然线性函数$f(\mathbf{x})=\alpha^T\mathbf{x}+\beta$ 既是凸函数又是凹函数，但都不严格。

证明很简单，但上述定理的逆定理不成立。我们称集合$H_S(f,c)$ 为函数$f$ 在集合$S$ 上关于$c$ 的水平集。

证明：

• 必要性。若$f$ 是凸函数，则有：

\begin{aligned} & \because f(\alpha\mathbf{y}+(1-\alpha)\mathbf{x})\leq \alpha f(\mathbf{y})+(1-\alpha)f(\mathbf{x})\\ & \therefore \frac{f(\alpha\mathbf{y}+(1-\alpha)\mathbf{x})-f(\mathbf{x})}{\alpha}\leq f(\mathbf{y})-f(\mathbf{x})\ \end{aligned}

然后把$f(\alpha\mathbf{y}+(1-\alpha)\mathbf{x})$ 在$\mathbf{x}$ 处泰勒展开：

$$f(\alpha\mathbf{y}+(1-\alpha)\mathbf{x})=f(\mathbf{x})+\nabla f(\mathbf{x})^T\alpha(\mathbf{y}-\mathbf{x})+\frac{\nabla^2 f(\zeta)}{2}\alpha^2(\mathbf{y}-\mathbf{x})^2$$

于是代入有：

$$\nabla f(\mathbf{x})^T(\mathbf{y}-\mathbf{x})+\frac{\nabla^2 f(\zeta)}{2}\alpha(\mathbf{y}-\mathbf{x})^2\leq f(\mathbf{y})-f(\mathbf{x})$$

令$\alpha\rightarrow 0$，则有$\nabla f(\mathbf{x})^T(\mathbf{y}-\mathbf{x})\leq f(\mathbf{y})-f(\mathbf{x})$。

• 充分性。假如已有上式，那么对于任意$\mathbf{x},\mathbf{y},\alpha$，我们取$\mathbf{z}=(1-\alpha)\mathbf{x}+\alpha\mathbf{y}$，代入上式：

$$\nabla f(\mathbf{z})^T(\mathbf{x}-\mathbf{z})\leq f(\mathbf{x})-f(\mathbf{z})\\ \nabla f(\mathbf{z})^T(\mathbf{y}-\mathbf{z})\leq f(\mathbf{y})-f(\mathbf{z})\\$$

把上两式加权和，有：

$$\alpha f(\mathbf{y})+(1-\alpha)f(\mathbf{x})\geq \nabla f(\mathbf{z})^T(\mathbf{z}-\mathbf{z})+f(\mathbf{z})=f(\mathbf{z})$$

证毕。

若一个 MP 问题的目标函数$f$ 是凸函数，可行集$X$ 是凸集，则称规划问题为凸规划。

证明：对于$g_i$ 是凸函数，根据前述定理，则它的零水平集 ($g_i(\mathbf{x})\leq 0$) 都为凸集。容易验证$h_j(\mathbf{x})=0$ 也是凸集，那么所有的凸集的交还是凸集。

证明：设$\mathbf{x}^*$ 是 MP 的一个局部最优解，故有一个邻域$N_\delta(\mathbf{x}^*)$ 满足：

$$\forall\mathbf{x}\in X\cap N_\delta(\mathbf{x}^*),f(\mathbf{x^*})\leq f(\mathbf{x})$$

反设$\mathbf{x}^*$ 不是最优解，那么存在一个$\bar{\mathbf{x}}\in X$，使得$f(\bar{\mathbf{x}})< f(\mathbf{x}^*)$，则有：

$$\begin{aligned} f(\alpha\bar{\mathbf{x}}+(1-\alpha)\mathbf{x}^*) & \leq \alpha f(\bar{\mathbf{x}})+(1-\alpha)f(\mathbf{x}^*)\\ &<\alpha f(\mathbf{x}^*)+(1-\alpha)f(\mathbf{x}^*)\\ & = f(\mathbf{x}^*) \end{aligned}$$

而当$\alpha$ 充分小时，$\alpha\bar{\mathbf{x}}+(1-\alpha)\mathbf{x}^*$ 可以成为邻域$N_\delta(\mathbf{x}^*)$ 中的点，故与$\mathbf{x}^*$ 是局部极小值相矛盾。证毕

# 一维搜索方法

考虑一个一维的单峰函数：$f(t),t\in[a,b]$，满足$\exists t^*$，使得$f$ 在$[a,t^*]$ 上单调增，在$[t^*,b]$ 上单调减。我们要解决的问题就是给出搜索区间$[a,b]$，去求得函数的峰$t^*$。

方法 1：三分法。
给定一个搜索区间$[a,b]$，我们可以在其中选择两个点$x_1,y_1$，满足$a<x_1<y_1<b$：

然后计算两个函数值$f(x_1),f(y_1)$。根据这两个函数值，我们一定可以保证峰不在两侧区间$[a,x_1]$ 或$[y_1,b]$ 中的某一个里。

很显然，若$f(x_1)>f(y_1)$，则说明函数的峰一定不在$[y_1,b]$ 中，反之若$f(x_1)<f(y_1)$ 则函数的峰比不可能在$[a,x_1]$ 中。这个结合图像特别好理解。因为若峰在$[a,x_1]$ 中，那么函数在$[x_1,b]$ 上就单调减，那么就不可能$f(x_1)<f(x_2)$，反之同理。

因此每分割一次区间，我们都可以去掉一部分，如此反复操作最后就可以把区间越变越小，逼近峰点。

---

方法 2：黄金分割法。

一个有点明显的优化就是我们不想每次都计算$f(x_1),f(y_1)$ 两个函数值，我们是否可以重复利用之前已经计算的结果？用图表示就是：

第一次搜索区间为$[a,b]$，我们选择$x_1,y_1$ 后扔掉了$[a,x_1]$ 区间。第二次搜索区间为$[x_1,b]$，我们选择$x_2,y_2$，此时我们希望$x_2=y_1$。这样就不用去计算$f(x_2)$，而只需利用之前计算的$f(y_1)$ 即可。

为保持这个过程可以永远进行下去，我们希望$x_1,y_1$ 在搜索区间$[a,b]$ 中的位置和比例关系和$y_1,y_2$ 在搜索区间$[x_1,b]$ 也$x,y$ 的选择是对称的（因为也可能扔掉右边的区间）：

故有比例关系：

$$1-\frac{y_1-a}{b-a}=\frac{y_1-x_1}{b-x_1}$$

设搜索区间$[a,b]$ 长度为单位 1，设$x_1-a=t$ 那么$b-y_2=t$，$y_1-x_1=1-2t$。那么代入就会有：

$$1-\frac{1-t}{1}=\frac{1-2t}{1-t}$$

又因为$0<t<1$ 解得$t=\frac{3-\sqrt{5}}{2}$。所以$y-a=\frac{\sqrt{5}-1}{2}(b-a)$。故分割过程可变为：

1. 选定$x_1=a+\frac{3-\sqrt{5}}{2}(b-a),y_1=a+\frac{\sqrt{5}-1}{2}(b-a)$ 计算$f(x_1),f(y_1)$。
2. 根据大小关系决定扔掉区间$[a,x_1]$ 或$[y_1,b]$。
3. 若扔掉$[a,x_1]$，则令搜索区间为$[x_1,b]$，且选定新的$x_2=y_1,y_2=x_1+\frac{\sqrt{5}-1}{2}(b-x_1)$。对称同理，如此反复。

---

方法 3： Newton 法

此时问题发生了变化，若$f$ 是一维可微的，即使它不是单峰的，我们也可以去尝试求解它的一个极值点。

假如我们现在已经有一个尝试点$t_k$（它不是极值点，但在前往极值点的路上 www），我们希望构建$t_{k+1}$，再这样一步步进行，直到逼近极值点。

首先，我们可以把原函数$f(x)$ 的值，都用在$t_k$ 处的二阶泰勒展开式近似：

$$\forall x,f(x)\approx f(t_k)+f'(t_k)(x-t_k)+\frac{f''(t_k)}{2}(x-t_k)^2$$

然后要求$f(x)$ 的极值点，就是让$f'(x)=0$，则有：

$$0=f'(x)\approx f'(t_k)+f''(t_k)(x-t_k)$$

解得：

$$x\approx t_k-\frac{f'(t_k)}{f''(t_k)}$$

当然这里是约等于！是因为显然由于我们只用了二阶近似，所以若$x=t_k-\frac{f'(t_k)}{f''(t_k)}$，那么$f'(x)\neq 0$。

但是这仍然有指导意义，因为$t_k-\frac{f'(t_k)}{f''(t_k)}$ 比$t_k$ 更接近极值点。故我们可以取新的尝试点$t_{k+1}=t_k-\frac{f'(t_k)}{f''(t_k)}$，如此反复计算以逼近极值点。

Newton 法显然不一定总是收敛到极值点的，实际上它具有局部收敛性，即尝试点若离极值点较近时会导致收敛。
在 https://subonan.com/2022/04/22 / 数值计算 / 中有给出了 Newton 法局部收敛性的定理说明。

# 无约束最优化方法

当可行域$X=R^n$ 时，我们讨论无约束的$n$ 元函数优化问题$min\;f(\mathbf{x})$(Unbounded Mathematical Programming, UMP )。

证明：泰勒展开，得

$$f(\bar{\mathbf{x}}+t\mathbf{p})=f(\bar{\mathbf{x}})+t\nabla f(\bar{\mathbf{x}})^T\mathbf{p}+o(\|t\mathbf{p}\|)$$

由于$\nabla f(\bar{\mathbf{x}})^T\mathbf{p}<0$，总可以取到足够小的$t$ 使得

$$t[\nabla f(\bar{\mathbf{x}})^T\mathbf{p}+\frac{o(\|t\mathbf{p}\|)}{t}]<0$$

故存在$t$，使得$f(\bar{\mathbf{x}}+t\mathbf{p})<f(\bar{\mathbf{x}})$。证毕

证明：当 Hessian 矩阵正定时，说明$f$ 是一个可微严格凸函数，所以有：

$$\forall \mathbf{x},\nabla f(\mathbf{x}^*)^T(\mathbf{x}-\mathbf{x}^*)< f(\mathbf{x})-f(\mathbf{x}^*)$$

又因为$\nabla f(\mathbf{x}^*)=0$，故有$\forall \mathbf{x},f(\mathbf{x}^*)<f(\mathbf{x})$。证毕

---

最速下降法

优化过程中，一个很显然的策略就是沿着梯度方向进行更新迭代。可是沿着哪个方向函数值下降最快呢？根据泰勒公式：

$$f(\bar{\mathbf{x}})-f(\bar{\mathbf{x}}+t\mathbf{p})=-t\nabla f(\bar{\mathbf{x}})^T\mathbf{p}-o(\|t\mathbf{p}\|)$$

假如$\nabla f(\bar{\mathbf{x}})$ 和$\mathbf{p}$ 长度一定的话，显然它们共线时，右式值最大，即

$$\mathbf{p}=-\nabla f(\bar{\mathbf{x}})$$

时，函数值下降最大。

---

共轭方向法

证明：假设存在一组不全为 0 的数满足：

$$\sum_{i=0}^{n-1}\alpha_i \mathbf{p}_i=0$$

用$\mathbf{p}_j^TA$ 左乘上式：

$$\sum_{i=0}^{n-1}\alpha_i \mathbf{p}_j^TA\mathbf{p}_i=0$$

由于$A$ 是正定的，且$\mathbf{p}_i$ 是共轭方向，故当且仅当$i=j$ 时，有：

$$\alpha_i \mathbf{p}_i^TA\mathbf{p}_i=0$$

其他情况$\mathbf{p}_j^TA\mathbf{p}_i=0$。又因为$A$ 正定，故$\mathbf{p}_i^TA\mathbf{p}_i>0$，故$\alpha_i=0$。同理可推出所有的$\alpha_i=0$。证毕。

考虑二次严格凸函数的无约束最优化问题 ( AP )：

$$min\;f(\mathbf{x})=\mathbf{x}^TA\mathbf{x}+\mathbf{b}^T\mathbf{x}+\mathbf{c}$$

其中$A$ 是$n$ 阶实对称正定矩阵。

对于 AP 问题，我们讨论形成一组共轭方向的通用方法：

1. 任选一个初始点$\mathbf{x}_0$，第一个搜索方向取$\mathbf{p}_0=-\nabla f(\mathbf{x}_0)$。
2. 沿着$\mathbf{p}_0$ 方向一维搜索，得到一个新的点：$\mathbf{x}_1=\mathbf{x}_0+t_0\mathbf{p}_0$。
3. 第二个搜索方向取$\mathbf{p}_1=-\nabla f(\mathbf{x}_1)+\lambda_0\mathbf{p}_0$。其中要求$\mathbf{p}_1$ 和$\mathbf{p}_0$ 相互共轭，利用$\mathbf{p}_1^T A\mathbf{p}_0=0$，解得：

$$\lambda_0=\frac{\mathbf{p}_0^TA\nabla f(\mathbf{x}_1)}{\mathbf{p}_0^TA\mathbf{p}_0}$$

4. 继续在$\mathbf{p}_1$ 方向上一维搜索。实际上，通式是：

$$\mathbf{p}_{k+1}=-\nabla f(\mathbf{x}_{k+1})+\lambda_k\mathbf{p}_k\\ \lambda_k=\frac{\mathbf{p}_k^TA\nabla f(\mathbf{x}_{k+1})}{\mathbf{p}_k^TA\mathbf{p}_k}$$

为了简化记忆，实际上：

$$\lambda_k=\frac{\|\nabla f(\mathbf{x}_{k+1})\|^2}{\|\nabla f(\mathbf{x}_k)\|^2}$$

以此产生了若干个搜索方向被称为 Fletcher-Reeves 共轭梯度法。

F-R 方法收敛较快，但极度依赖于一维搜索的精度。当一维搜索并不保证精度时，我们可以用 Polak-Ribiere-Polyak 共轭梯度法，区别仅仅在于：

$$\lambda_k=\frac{\nabla f(\mathbf{x}_{k+1})^T[\nabla f(\mathbf{x}_{k+1})-\nabla f(\mathbf{x}_k)]}{\|\nabla f(\mathbf{x}_k)\|^2}$$

# 约束最优化方法

# K-T 条件法

设 MP 问题得可行域为：

$$X=\{\mathbf{x}\in R^n\;|\;g_i(\mathbf{x})\leq 0,h_j(\mathbf{x})=0\},i=1,...,p,j=1,...,q$$

对于$g_i(\mathbf{x})\leq 0$ 的约束也可以继续分类，我们称对于使得$g_i(\mathbf{x}^*)=0$ 的点$\mathbf{x}^*$，$g_i(\mathbf{x}^*)=0$ 是它的一个积极约束，一个微小的扰动就会导致积极约束被破坏。我们记$\mathbf{x}^*$ 的积极约束的下标为：

$$I(\mathbf{x}^*)=\{i\;|\;g_i(\mathbf{x}^*)=0\}$$

其中向量组$\nabla g_i(\mathbf{x}^*),i\in I(\mathbf{x}^*),\nabla h_j(\mathbf{x}^*)$ 线性无关称为一个约束规范条件。有各种各样的约束规范条件，当然这个是比较简单的一种。

式子$\nabla f(\mathbf{x}^*)+\sum_{i\in I(\mathbf{x}^*)}\lambda_i^*\nabla g_i(\mathbf{x}^*)+\sum_j\mu_j^*\nabla h_j(\mathbf{x}^*)=0$ 被称为 MP 的 Kuhn-Tucker 条件，凡是满足 K-T 条件的点被称为 MP 的 K-T 点。上述定理说明， MP 的局部最优解一定是 MP 的 K-T 点。

下面我们来考虑一下 K-T 条件的几何意义。首先若不考虑$h_j$ 的部分，则有：

$$-\nabla f(\mathbf{x}^*)=\sum_{i\in I(\mathbf{x}^*)}\lambda_i^*\nabla g_i(\mathbf{x}^*),\quad \lambda_i^*\ge 0$$

可以发现负梯度$-\nabla f(\mathbf{x}^*)$ 落在了由$\nabla g_i(\mathbf{x}^*),i\in I(\mathbf{x}^*)$ 的线性组合的集合中：

$$C=\{\mathbf{y}\;|\;\mathbf{y}=\sum_{i\in I(\mathbf{x}^*)}\alpha_i\nabla g_i(\mathbf{x}^*),\alpha_i\geq 0\}$$

上述集合称为由$\nabla g_i(\mathbf{x}^*)$ 张成的锥。注意，由于$\alpha_i\geq 0$，所以是锥而不是子空间。

相比之下，若考虑$h_j$ 的部分，则有：

$$-\nabla f(\mathbf{x}^*)=\sum_j\mu_j^*\nabla h_j(\mathbf{x}^*)$$

这只能说明负梯度落在子空间内，而不是锥内。此外，这个条件就是拉格朗日乘子法中的条件。

若在 K-T 条件中，要求$g_i(\mathbf{x}),i\in I$ 在$\mathbf{x}^*$ 处均可微，那么 K-T 条件可进一步为：

$$\nabla f(\mathbf{x}^*)+\sum_{i=1}^p\lambda_i^*\nabla g_i(\mathbf{x}^*)+\sum_j\mu_j^*\nabla h_j(\mathbf{x}^*)=0\\ \lambda_i^*\geq 0,i=1,...,p\\ \lambda_i^*g_i(\mathbf{x}^*)=0,i=1,...,p$$

为什么现在不需要区分$I(\mathbf{x}^*)$ 了呢？就是因为有了互补松紧条件 ($\lambda_i^* g_i(\mathbf{x}^*)=0$)，若$i\notin I(\mathbf{x}^*$)，则$\lambda_i^*=0$。

引入 Lagrange 函数：

$$L(\mathbf{x},\lambda,\mu)=f(\mathbf{x})+\sum_{i=1}^p\lambda_ig_i(\mathbf{x})+\sum_{j=1}^q\mu_jh_j(\mathbf{x})$$

MP 问题的 K-T 条件可写为：

$$\begin{cases} \nabla_x L(\mathbf{x}^*,\lambda^*,\mu^*)=0\\ \lambda_i^*g_i(\mathbf{x}^*)=0,i=1,..,p\\ \lambda_i^*\geq 0,i=1,...,p \end{cases}$$

其中，$\lambda,\mu$ 称为 Lagrange 乘子。

所以一个基本的求解方法就可以为：

1. 确定$f,g,h$ 连续可微性（实际中一般用的函数都是处处连续可微的）。
2. 确定$f,g$ 是否凸，$h$ 是否线性。
3. 列出 K-T 条件和可行性条件，并进行求解。
4. 检验解是否满足上述定理的要求。

---

# 简约梯度法

考虑如下问题：

$$\begin{cases} min & f(\mathbf{x})=0\\ s,t. & A\mathbf{x}=\mathbf{b}\\ & \mathbf{x}\geq 0 \end{cases}$$

其中$A$ 是一个秩为$m$ 的$m\times n$ 的矩阵。现作如下约束非退化假设：

1. 每个可行点至少有$m$ 个大于等于 0 的分量。
2. $A$ 的任意$m$ 列线性无关。

类似之前的单纯形法，我们仍可以把$\mathbf{x}$ 拆分为基向量和非基向量两部分，即：

$$f(\mathbf{x})=f(\mathbf{x}_B,\mathbf{x}_N)$$

仍有$B\mathbf{x}_B+N\mathbf{x}_N=\mathbf{b}$，故$f$ 可以表示为关于$\mathbf{x}_N$ 的函数：

$$F(\mathbf{x}_N)=f(\mathbf{x}_B,\mathbf{x}_N)=f(B^{-1}\mathbf{b}-B^{-1}N\mathbf{x}_N,\mathbf{x}_N)$$

对它求梯度，有：

$$\nabla F(\mathbf{x}_N)=-(B^{-1}N)^T\nabla_B f(\mathbf{x})+\nabla_N f(\mathbf{x})$$

上述梯度称为$f$ 在$\mathbf{x}$ 处的简约梯度。其中，$\nabla_B f(\mathbf{x})$ 是$f$ 相对于那些基分量求偏导组成的向量，$\nabla_N f(\mathbf{x})$ 是对非基分量求偏导组成的向量。即有$\nabla f(\mathbf{x})=\begin{pmatrix}\nabla_B f(\mathbf{x})\\ \nabla_N f(\mathbf{x})\end{pmatrix}$。

---

下面我们考虑如何构造搜索方向$\mathbf{p}=\begin{pmatrix}\mathbf{p}_B\\ \mathbf{p}_N\end{pmatrix}$, 对于一个$\mathbf{x}$ 的最大的$m$ 个分量组成的向量为$\mathbf{x}_B>0$（非退化假设），并记下标集合为$I_B$。

一个朴素的想法就是$\mathbf{p}_N=-\nabla F(\mathbf{x}_N)$。这样作会导致$\mathbf{x}\geq 0$ 的条件被破坏。想象一下如果沿着$-\nabla F(\mathbf{x}_N)$ 的方向一维搜索，那$\mathbf{x}_N$ 中如果有一个分量为 0，那么对于$\mathbf{x}_N-t\nabla F(\mathbf{x}_N)$，我们就无法取足够小的$t$ 使其为正。

我们这样选取$\mathbf{p}_N$ 的分量：

$$\forall i\notin I_B\\ p_i=\begin{cases} -\nabla F(\mathbf{x}_N)_i & \nabla F(\mathbf{x}_N)_i\leq 0\\ -x_i\nabla F(\mathbf{x}_N)_i & \nabla F(\mathbf{x}_N)_i>0 \end{cases}$$

这样可以巧妙地避免上述问题，因为$x_i(1-t\nabla F(\mathbf{x}_N)_i)$ 的话，若$x_i=0$ 则无论怎么取$t$ 也不影响，若$x_i>0$，则可以选择足够小的$t$ 使得$1-t\nabla F(\mathbf{x}_N)_i>0$。

当然另一方面我们还希望保持$A\mathbf{x}=\mathbf{b}$ 的可行性，即：

$$\forall t,A(\mathbf{x}+t\mathbf{p})=\mathbf{b}$$

那么我们就想让$A\mathbf{p}=0$。就可以解得：

$$\mathbf{p}_B=-B^{-1}N\mathbf{p}_N$$

如此我们可以定下，若给了尝试点$\mathbf{x}$ 满足$A\mathbf{x}=\mathbf{b},\mathbf{x}\geq 0$，构造一个搜索方向为：

$$\begin{cases} \mathbf{p}_N=\begin{cases} -\nabla F(\mathbf{x}_N)_i & \nabla F(\mathbf{x}_N)_i\leq 0\\ -x_i\nabla F(\mathbf{x}_N)_i & \nabla F(\mathbf{x}_N)_i>0 \end{cases}\\ \mathbf{p}_B=-B^{-1}N\mathbf{p}_N \end{cases}$$

证明：
关于 1，可行性在构造过程中已经说明了。关于下降，其实有：

$$\nabla f(\mathbf{x})^T\mathbf{p}=\nabla F(\mathbf{x}_N)^T\mathbf{p}_N\leq 0$$

由本章 "无约束优化" 的第一个定理，它是是一个函数下降方向。

关于 2，我们先证明若点$\mathbf{x}$ 满足 K-T 条件，则$\mathbf{p}=0$。这个一次的$A\mathbf{x}=\mathbf{b}$ 满足 K-T 条件，即：

$$\begin{cases} \nabla f(\mathbf{x})-\lambda^*+A^T\mu^*=0\\ \lambda^{*T}\mathbf{x}=0\\ \lambda^*\geq 0 \end{cases}$$

同时我们可以把$\lambda^*$ 分解为$(\lambda_B^*,\lambda_N^*)$，然后因为$\mathbf{x}_B>0$，由互补松紧性，有$\lambda_B^*=0$。又因为：

$$\nabla_B f(\mathbf{x})-\lambda_B^*+B^T\mu^*=0$$

注意！这里逻辑在于并不是$A+B=0\Rightarrow A=0\wedge B=0$，而是若一个向量为 0，则它的所有分量都为 0。所以我们可以从$\nabla f(\mathbf{x})-\lambda^*+A^T\mu^*=0$ 推出上面这个式子。

再根据简约梯度的式子，有：

$$\begin{aligned} &\because\nabla f(\mathbf{x})-\lambda^*+A^T\mu^*=0\\ &\therefore\nabla_B f(\mathbf{x})-\lambda_B^*+B^T\mu^*=0\\ &\because \lambda_B^*=0\\ &\therefore\mu^*=-(B^T)^{-1}\nabla_B f(\mathbf{x})\\ &\because\nabla f(\mathbf{x})-\lambda^*+A^T\mu^*=0\\ &\therefore\nabla_N f(\mathbf{x})-\lambda_N^*+N^T\mu^*=0\\ &\therefore \lambda_N^*=\nabla F(\mathbf{x}_N)\\ \end{aligned}$$

又因为$\lambda^{*T}\mathbf{x}=0,\lambda^*\geq 0$，根据$\mathbf{p}$ 的构造方法，故有$\mathbf{p}_N=0$。所以$\mathbf{p}_B$ 也为 0。

反过来的话，很显然我们可以取$\lambda_B^*=0,\lambda_N^*=\nabla F(\mathbf{x}_N),\mu^*=-(B^T)^{-1}\nabla_B f(\mathbf{x})$ 即可验证满足 K-T 条件。

证毕

---

有了以上搜索方向，我们再回看一下沿着这个方向一维搜索的过程。给定点$\mathbf{x}$ 和搜索方向$\mathbf{p}$，我们希望:

$$\forall t>0,\mathbf{x}+t\mathbf{p}\geq 0$$

于是$t$ 的上界为：

$$\begin{cases} +\infty & \mathbf{p}\geq 0\\ \mathop{min}\limits_{1\leq i\leq n}\{-\frac{x_i}{p_i}\;|\;p_i<0\}&otherwise \end{cases}$$

---

在非退化假设下，上述 Wolfe 简约梯度方法是收敛的。我们总结它的计算步骤为：

1. 选取初始可行点$\mathbf{x}$。
2. 设$I_B$ 是$\mathbf{x}$ 的最大的$m$ 个分量的下标集合，$A=(B,N)$。计算：

$$\nabla f(\mathbf{x})=\begin{pmatrix} \nabla_B f(\mathbf{x})\\ \nabla_N f(\mathbf{x}) \end{pmatrix}$$

3. 计算简约梯度：

$$\nabla F(\mathbf{x}_N)=-(B^{-1}N)^T\nabla_B f(\mathbf{x})+\nabla_N f(\mathbf{x})$$

4. 用下述方法构造搜索方向$\mathbf{p}$：

$$\begin{cases} \mathbf{p}_N=\begin{cases} -\nabla F(\mathbf{x}_N)_i & \nabla F(\mathbf{x}_N)_i\leq 0\\ -x_i\nabla F(\mathbf{x}_N)_i & \nabla F(\mathbf{x}_N)_i>0 \end{cases}\\ \mathbf{p}_B=-B^{-1}N\mathbf{p}_N \end{cases}$$

5. 进行一维搜索，确定$min\;f(\mathbf{x}+t\mathbf{p})$。
6. 重复 2-5 直到满足停止条件。

# 惩罚函数法

一个自然的想法就是，对于一个带约束的优化问题：$\mathop{min}\limits_{\mathbf{x}\in X}f(\mathbf{x})$，我们可以设置一个罚函数：

$$p_c(\mathbf{x})=\begin{cases} 0 & \mathbf{x}\in X\\ c & otherwise \end{cases}$$

其中$c$ 反映了惩罚力度，然后把原问题转化为无约束的规划问题$min f(\mathbf{x})+p_c(\mathbf{x})$。

但是这样不好的是显然这个罚函数有跳跃，即并不连续可微。取而代之地可以选择一个新的罚函数：

$$p_c(\mathbf{c})=c\sum_{i=1}^pmax(g_i(\mathbf{x}),0)^2+\frac{c}{2}\sum_{j=1}^qh_j(\mathbf{x})^2$$

可以证明，当$f,g_i,h_i$ 均连续可微时，上述罚函数也是连续可微的。

然而在实际过程中，一下就选取到合适的惩罚因子$c$ 是很困难的。相反，我们会构造一个等比数列：$\{c_0,\alpha c_0,\alpha^2c_0,...\}$ 这样同时求解若干个无约束优化问题。又因为乘法因子等比，所以在并行地求解过程中，乘法函数也很好计算$p_{\alpha c_0}=\alpha p_{c_0}$，即$\{p_{c_0}(\mathbf{x}),\alpha p_{c_0}(\mathbf{x}),\alpha^2p_{c_0}(\mathbf{x}),...\}$ 可一起计算。最后在所有惩罚因子中选取最小的那个。

此外，我们也可以选择障碍函数。仅考虑含不等式的约束$g_i(\mathbf{x})\leq 0$ 的话，我们可以设置障碍函数：

$$B_c(\mathbf{x})=-c\sum_{i=1}^p\frac{1}{g_i(\mathbf{x})}$$

故趋近于边界时，总会有一个$g_i(\mathbf{x})$ 从负数趋近于 0，即障碍函数会趋近于无穷大。

由于障碍函数增长速度很快，在实际中我们往往要构造一个单调减小的等比数列，然后求得极优值。