你怎么不牵我!

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

线性规划问题的规范形式

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

线性规划问题的标准形式

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

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

我们考虑标准形式的 LP,假设rank(A)=mrank(A)=mAAm×nm\times n 的),故AA 中有mm 个线性无关的列向量,它们构成了满秩方阵Bm×mB_{m\times m}。再把AA 中其余各列
组成的子阵记为NN,即A=(B,N)A=(B,N),再把x=(x1,x2,...,xn)T\textbf{x}=(x_1,x_2,...,x_n)^T 分为两部分:记为xB\textbf{x}_BxN\textbf{x}_N,则有:

BxB+NxN=bB\textbf{x}_B+N\textbf{x}_N=\textbf{b}

可以解得

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

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

定义:设BB 是秩为mm 的约束矩阵ARm×nA\in R^{m\times n} 中的一个mm 阶满秩方阵,则称BB 为一个基(或基阵)BBmm 个线性无关的列向量称为基向量,变量x\textbf{x} 中与之对应的mm 个分量称为基变量,其余分量称为非基变量。令所有非基变量为 0,得到的解x=(xBxN)=(B1b0)\textbf{x}=\begin{pmatrix}\textbf{x}_B\\ \textbf{x}_N\end{pmatrix}=\begin{pmatrix}B^{-1}\textbf{b}\\0\end{pmatrix} 称为相对BB基本解,当B1b0B^{-1}\textbf{b}\geq 0 时,称基本解x\textbf{x}基本可行解,这时对应的基BB 称为可行基

定理:可行解x\textbf{x} 是基本可行解,当且仅当它的正的分量对应的AA 中的列向量线性无关。

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

反之,由于x\textbf{x} 是可行解,故有x0\textbf{x}\geq 0Ax=bA\textbf{x}=\textbf{b}。不失一般性假设x\textbf{x} 的前kk 个分量为正:

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

那么显然有kmk\leq m,可以讨论:

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

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

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

定理:一个标准形式的 LP 问题若有可行解,则至少有一个基本可行解。

证:设x0=(x10,...,xn0)T\textbf{x}^0=(x_1^0,...,x_n^0)^T 是一个可行解,则有Ax0=b,x00A\textbf{x}_0=\textbf{b},\textbf{x}^0\geq 0。不失一般性,设x0\textbf{x}^0 的正分量为前kk 个。如果它们对应的AA 中的列向量A1,...,AkA_1,...,A_k 已经线性无关了,则根据前述定理,x0\textbf{x}^0 已经是一个基本可行解了。

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

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

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

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

由于xj0>0,j=1,...,kx_j^0>0,j=1,...,k,故可以选择一个足够小的正数ε\varepsilon 使得:

x0±εδ0\textbf{x}^0\pm\varepsilon\delta\geq 0

故有

A(x0±εδ)=Ax0+εAδ=bA(\textbf{x}^0\pm\varepsilon\delta)=A\textbf{x}^0+\varepsilon A\delta=\textbf{b}

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

定理:一个标准形式的 LP 问题若有有限的最优值,则一定存在一个基本可行解是最优解。

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

cT(x0εδ)=cTx0εcTδcT(x0+εδ)=cTx0+εcTδ\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

因为cTx0\textbf{c}^T\textbf{x}^0 是最优解,故有:

cT(x0±εδ)cTx0\textbf{c}^T(\textbf{x}^0\pm\varepsilon\delta)\geq\textbf{c}^T\textbf{x}^0\\

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

# LP 问题的单纯形法

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

假设已经找到了一个非退化的基本可行解xˉ\bar{\textbf{x}},即找到了一个基BB,此时可以将方程组Ax=bA\textbf{x}=\textbf{b} 转化为与之等价的方程组:

xB+B1NxN=B1b\textbf{x}_B+B^{-1}N\textbf{x}_N=B^{-1}\textbf{b}

不失一般性假定B=(A1,...,Am)B=(A_1,...,A_m),记向量:

Ajˉ=B1Ajbˉ=B1b\bar{A_j}=B^{-1}A_j\\ \bar{\textbf{b}}=B^{-1}\textbf{b}

则可以把方程组转化为:

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

这个方程组被称为对应于基BB典则方程组,简称典式
相应的,可以改写目标函数值为:

cTx=cBTxB+cNTxN=cBT(bˉB1NxN)+cNTxN=cBTbˉj=m+1n(cBTAˉjcj)xj\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}

其中不失一般性地假设了x\textbf{x} 的前mm 个分量对应xB\textbf{x}_B,后nmn-m 个分量对应了xN\textbf{x}_N。可以注意到:

(Aˉ1,...,Aˉm)=B1(A1,...,Am)=B1B=I(\bar{A}_1,...,\bar{A}_m)=B^{-1}(A_1,...,A_m)=B^{-1}B=I

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

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

引入记号:

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

向量表示为:

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

从而目标函数值为:

cTx=cBTbˉζTx\textbf{c}^T\textbf{x}=\textbf{c}_B^T\bar{\textbf{b}}-\zeta^T\textbf{x}

定理 (最优性准则):如果ζ0\zeta\leq 0,则xˉ=(bˉ0)\bar{\textbf{x}}=\begin{pmatrix}\bar{\textbf{b}}\\0\end{pmatrix} 是原问题的最优解。

证:对于任意一个可行解x\textbf{x},一定有x0\textbf{x}\geq 0,那么就有:

cTxˉ=cBTbˉcBTbˉζTx=cTx\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}

定理:如果向量ζ\zeta 的第kk 个分量ζk>0\zeta_k>0(显然m+1knm+1\leq k\leq n),且向量Aˉk=B1Ak0\bar{A}_k=B^{-1}A_k\leq 0,则原问题无界。

证:令d=(Aˉk0)+ek\textbf{d}=\begin{pmatrix}-\bar{A}_k\\0\end{pmatrix}+e_k,其中eke_k 是第kk 个分量为 1,其余分量为 0 的列向量。
因为Aˉk0\bar{A}_k\leq 0,所以有d0\textbf{d}\geq 0。而:

Ad=(B,N)(Aˉk0)+Aek=Ak+Ak=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,观察xˉ+θd\bar{\textbf{x}}+\theta\textbf{d},此时有:

A(xˉ+θd)=Axˉ+θAd=bxˉ+θd0A(\bar{\textbf{x}}+\theta\textbf{d})=A\bar{\textbf{x}}+\theta A\textbf{d}=\textbf{b}\\ \bar{\textbf{x}}+\theta\textbf{d}\geq 0

所以xˉ+θd\bar{\textbf{x}}+\theta\textbf{d} 是原问题的可行解,那么目标函数值为:

cT(xˉ+θd)=cTxˉ+θ(cBT,cNT)(Aˉk0)+θcTek=cTxˉθ(cBTAˉkck)=cTxˉθζk\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}

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

定理:对于非退化的基本可行解xˉ\bar{\textbf{x}},若有ζ\zeta 中有任一分量ζk>0\zeta_k>0,且其对应的列向量Aˉk\bar{A}_k 至少有一个正分量,则能找到一个新的基本可行解x^\hat{\textbf{x}},使得cTx^<cTxˉ\textbf{c}^T\hat{\textbf{x}}<\textbf{c}^T\bar{\textbf{x}}

证:只需要找出x^\hat{\textbf{x}} 是什么。
类似前述证明,令

d=(Aˉk0)+ek\textbf{d}=\begin{pmatrix}-\bar{A}_k\\0\end{pmatrix}+e_k

则有Ad=0A\textbf{d}=\textbf{0}。令:

x^=xˉ+θd=(bˉθAˉk0)+θek\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。首先为了x^0\hat{\textbf{x}}\geq 0, 可以选

θ=min{bˉiaˉikaˉik>0,i=1,...,m}=bˉraˉrk\theta=min\{\frac{\bar{b}_i}{\bar{a}_{ik}}|\bar{a}_{ik}>0,i=1,...,m\}\\ =\frac{\bar{b}_r}{\bar{a}_{rk}}

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

x^i=bˉibˉraˉrkaˉik,i=1,...,m,irx^r=0x^k=θ=bˉraˉrkx^j=0,j=1,...,n,jk\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}

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

A1,...,Ar1,Ak,Ar+1,...,AmA_1,...,A_{r-1},A_k,A_{r+1},...,A_m

是线性无关的。(即把第rr 个分量换出去,把第kk 个换进来)。
反设不是线性无关的,由于先前的A1,...,Ar,...,AmA_1,...,A_r,...,A_m 是线性无关的,那么就会有:

Ak=i=1,irmyiAiA_k=\sum_{i=1,i\neq r}^my_iA_i

同时有:

Ak=BAˉk=i=1maˉikAiA_k=B\bar{A}_k=\sum_{i=1}^m\bar{a}_{ik}A_i

将上述两式相减可得:

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

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

cTx^=cTxˉθζk<cTxˉ\textbf{c}^T\hat{\textbf{x}}=\textbf{c}^T\bar{\textbf{x}}-\theta\zeta_k<\textbf{c}^T\bar{\textbf{x}}

证毕。

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

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

定理:对于任何一个非退化的线性规划问题,从任何基本解开始,经过有限次换基迭代,或得到一个最优的基本可行解,或作出该问题无界的判定。

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

单纯形法步骤:

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

# LP 问题的初始解

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

{min  z=cTxs.t.Ax=b(b0)x0\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,b)(A,\textbf{b}) 行变换,然后行变换出来一个单位矩阵后,就得到了一个基本可行解。但实际上无法确认究竟行变换后b\textbf{b} 是否能保证不为负数。譬如我不能先入为主就认为x1,x2,x30,x4,x5=0x_1,x_2,x_3\neq 0,x_4,x_5=0 是一组可行解,很可能做完行变换后b\textbf{b} 就不是非负的了。

为了解决 "得到一个初始的基本可行解" 这个问题,我们引入mm 个人工变量:

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

并研究以下 LP 问题:

{min  g=i=n+1n+mxis.t.Ax+xα=bx0,xα0\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=0g=0。此时即xα=0\textbf{x}_\alpha=0
而上述 LP 问题已经有了一个天然的初始基本可行解,即后mm 个变量。

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

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

# LP 问题的对偶性

定义:一个规范形式的 LP 问题:

{min  cTxs.t.Axbx0\begin{cases} min\; & \textbf{c}^T\textbf{x}\\ s.t. & A\textbf{x}\geq \textbf{b}\\ & \textbf{x}\geq 0 \end{cases}

它的对偶问题为:

{max  bTys.t.ATycy0\begin{cases} max\; & \textbf{b}^T\textbf{y}\\ s.t. & A^T\textbf{y}\leq \textbf{c}\\ & \textbf{y}\geq 0 \end{cases}

当一个 LP 问题为标准形式:

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

它的对偶问题为:

{max  bTys.t.ATyc\begin{cases} max\; & \textbf{b}^T\textbf{y}\\ s.t. & A^T\textbf{y}\leq \textbf{c}\\ \end{cases}

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

{min  2x1+3x2s.t.x1+x21x22x1,x20\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 加上第二个不等式就得到了2x1+3x242x_1+3x_2\geq 4
我们就得到了一个目标函数的一个下界,而这个下界是可能取不到的,譬如这个例子就是。
但是目标函数一定是大于等于这个下界的。
所以对偶问题实际上就是在寻找一个原问题的下界。

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

定理:如果一个 LP 问题有最优解,那么它的对偶问题也有最优解,且它们最优解相等。

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

Ax=byTAx=yTbATyc,  x0yTAxcTxyTb=yTAxcTx\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}

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

ATyc=AT(cB^TB^1)Tc=ζ0A^T\textbf{y}-\textbf{c}=A^T(\textbf{c}_{\hat{B}}^T\hat{B}^{-1})^T-\textbf{c}=\zeta\leq 0\\

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

yTb=cB^TB^1b=cB^Tx^B\textbf{y}^T\textbf{b}= \textbf{c}_{\hat{B}}^T\hat{B}^{-1}\textbf{b}=\textbf{c}_{\hat{B}}^T\hat{\textbf{x}}_B

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

定理:若x,w\textbf{x},\textbf{w} 分别是原问题及其对偶问题的可行解,那么它们分别是原问题、对偶问题的最优解的充要条件为cTx=wTb\textbf{c}^T\textbf{x}=\textbf{w}^T\textbf{b}

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

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

定理:对偶问题的对偶问题就是原问题。

定理 (互补松紧性):设x,w\textbf{x},\textbf{w} 分别是原问题和对偶问题的可行解,那么它们分别是原问题、对偶问题的最优解的充要条件是:

yT(Axb)=0xT(cATy)=0 \textbf{y}^T(A\textbf{x}-\textbf{b})=0\\ \textbf{x}^T(\textbf{c}-A^T\textbf{y})=0

证:考虑规范形式,就会有:

x,y,(Axb),(cATy)0yT(Axb)0,xT(cATy)0yT(Axb)+xT(cATy)=xTcyTb\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}

所以xTc=yTbyT(Axb)=0,xT(cATy)=0\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算法 却是从多面体内部一个点开始,产生一个直接穿过多面体内部的点列:
    1
    Karmarkar 证明了该算法迭代次数的上届是O(nL)O(nL) 每次迭代的平均计算量为O(n2.5)O(n^{2.5}),因而其算法的总的时间复杂度为O(n3.5L)O(n^{3.5}L),其中nn 是维数,LL 为输入规模。
  • 之后也诞生了诸如仿射尺度法、路径追踪法等 LP 问题的解法,其中最好的时间复杂度在O(n3L)O(n^3L)。已经在商业化软件中实现,但是内点法有一个固有弱点就是仅提供近似最优解。需另外加一个纯化的过程。
  • 单纯形法有 “热启动” 的特性,即从最后求得的基再开始能大大缩短求解时间。

# 整数线性规划

# 整数线性规划问题 ILP

定义整数规划问题 (Integer Linear Programming, ILP ) 是指下述规划问题:

{min  z=cTxs.t.  Ax=bx0xiI,iJ{1,2,,n}\begin{cases} min\;& z=\mathbf{c}^T\mathbf{x}\\ s.t.\;& A\textbf{x}=\textbf{b}\\ & \textbf{x}\geq 0\\ & x_i\in I, i\in J\subset\{1,2,\cdots,n\} \end{cases}

其中,AAm×nm\times n 的矩阵,cRn,bRm,x=(x1,...,xn)T,I\mathbf{c}\in R^n,\mathbf{b}\in R^m,\mathbf{x}=(x_1,...,x_n)^T,I 是自然数集合。

当我们想求解一个 ILP 问题时,可以想象先求解对应的 LP 问题,然后自然地把结果舍入到一个整数解。但是从 LP 问题的可行解舍入到一个可行的整数解是困难的,甚至是不可行的:
2
实际上这个舍入的过程复杂度已经相当于求解原问题。

# Gomory割平面法

考虑一种特殊的整数线性规划问题,即纯整数线性规划问题 (PP):

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

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

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

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

  1. DD0D\subset D_0
  2. P0P_0 无可行解,则PP 无可行解。
  3. P0P_0 的最优值是PP 最优值的一个下界。
  4. P0P_0 的最优解是一个整数向量,则他也是PP 的最优解。

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

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


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

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

xBi+jSˉaˉijxj=biˉ,i=1,...,mx_{B_i}+\sum_{j\in\bar{S}}\bar{a}_{ij}x_j=\bar{b_i},i=1,...,m

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

由于上述典式中,xBi,xjx_{B_i},x_j 都是非负的,因此有:

jSˉ[aˉij]xjjSˉaˉijxjxBi+jSˉ[aˉij]xjbiˉ,i=1,...,m\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 问题的可行域不变。很显然对于上述下取整后的式子,由于左侧都是整数,那么我们对右侧也下取整,并不影响不等式:

xBi+jSˉ[aˉij]xj[bˉi],i=1,...,mx_{B_i}+\sum_{j\in\bar{S}}[\bar{a}_{ij}]x_j\leq[\bar{b}_i],i=1,...,m

将原来的等式和上述不等式相减,得到:

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

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

jSˉfljxj+s=fl-\sum_{j\in\bar{S}}f_{lj}x_j+s=-f_l

定理:增加上述条件后,原 ILP 问题的可行域不变,但是x0\mathbf{x}^0 将不再可行。

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

xT=(xB1,...,xBm,s,xSˉ1,...,xSˉn)=(bˉ1,...,bˉm,fl,0,...,0)\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}

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


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

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

# 分枝定界法

还是考虑纯整数规划问题 (PP):

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

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

{min  z=cTxs.t.  Ax=bx0xl[xl]{min  z=cTxs.t.  Ax=bx0xl[xl]+1\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}

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

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

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

# 非线性规划

# 基本概念

定义:令xT=(x1,...,xn)\mathbf{x}^T=(x_1,...,x_n)nn 维欧氏空间RnR^n 中的一个点。f(x),gi(x),i=1,...,p,hj(x),j=1,...,qf(\mathbf{x}),g_i(\mathbf{x}),i=1,...,p,h_j(\mathbf{x}),j=1,...,q 是定义在RnR^n 上的实值函数。我们称如下形式的数学模型为数学规划(Mathematical Programming, MP ):

{minf(x)s.t.gi(x)0,i=1,...,phj(x)=0,j=1,...,q\begin{cases} min & f(\mathbf{x})\\ s.t. & g_i(\mathbf{x})\leq 0,i=1,...,p\\ & h_j(\mathbf{x})=0,j=1,...,q\\ \end{cases}

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

称为 MP可行域约束集。故 MP 问题可简记为minxXf(x)\mathop{min}\limits_{\mathbf{x}\in X}f(\mathbf{x})。当f,gi,hjf,g_i,h_j 函数中为非线性函数时,就将原 MP 称为非线性规划问题。

定义:对于非线性规划 MP 问题,若xX\mathbf{x}^*\in X,并且存在x\mathbf{x}^* 的一个邻域Nδ(x)={xRn    xx<δ}N_\delta(\mathbf{x}^*)=\{x\in R^n\;|\;\|\mathbf{x}-\mathbf{x}^*\|<\delta\},使得:

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

则称x\mathbf{x}^*MP 问题的局部最优点局部最小点,称f(x)f(\mathbf{x}^*) 称为 MP 问题的局部最优值局部极小值

若有f(x)<f(x)f(\mathbf{x}^*)<f(\mathbf{x}),那么称x\mathbf{x}^*MP 问题的严格局部最优点严格局部最小点,称f(x)f(\mathbf{x}^*) 称为 MP 问题的严格局部最优值严格局部极小值

# 凸函数和凸规划

定义:设SRnS\subset R^nnn 维欧氏空间中的一个点集,若对任意x,yS,λ[0,1]\mathbf{x},\mathbf{y}\in S,\lambda\in[0,1] 都有λx+(1λ)yS\lambda\mathbf{x}+(1-\lambda)\mathbf{y}\in S,就称SS凸集

若有实值函数f:SR1f:S\rightarrow R^1 且对任意的α(0,1)\alpha\in(0,1) 都有:

x,yS,f(αx+(1α)y)αf(x)+(1α)f(y)\forall \mathbf{x},\mathbf{y}\in S,f(\alpha\mathbf{x}+(1-\alpha)\mathbf{y})\leq \alpha f(\mathbf{x})+(1-\alpha)f(\mathbf{y})

则称ffSS 上的凸函数。若满足严格小于而非小于等于,则称为严格凸函数。若f-fSS 上的 (严格) 凸函数,则称ffSS 上的 (严格) 凹函数

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

定理:凸函数的数乘αf\alpha f、加和f1+f2f_1+f_2 都是凸函数,但凸函数的乘积不一定是凸函数。

定理:设ff 是凸函数,则集合:

HS(f,c)={xS    f(x)c}H_S(f,c)=\{\mathbf{x}\in S\;|\;f(\mathbf{x})\leq c\}

是凸集。

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

定理:设SRnS\subset R^n 是非空开凸集,f:SR1f:S\rightarrow R^1 可微,则:

ffSS 上的凸函数的充要条件是:

x,yS,f(x)T(yx)f(y)f(x)\forall \mathbf{x},\mathbf{y}\in S,\nabla f(\mathbf{x})^T(\mathbf{y}-\mathbf{x})\leq f(\mathbf{y})-f(\mathbf{x})

若要求ff 严格凸,那么充要条件需改为严格小于。

证明:

  • 必要性。若ff 是凸函数,则有:
\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(αy+(1α)x)f(\alpha\mathbf{y}+(1-\alpha)\mathbf{x})x\mathbf{x} 处泰勒展开:

f(αy+(1α)x)=f(x)+f(x)Tα(yx)+2f(ζ)2α2(yx)2f(\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

于是代入有:

f(x)T(yx)+2f(ζ)2α(yx)2f(y)f(x)\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})

α0\alpha\rightarrow 0,则有f(x)T(yx)f(y)f(x)\nabla f(\mathbf{x})^T(\mathbf{y}-\mathbf{x})\leq f(\mathbf{y})-f(\mathbf{x})

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

f(z)T(xz)f(x)f(z)f(z)T(yz)f(y)f(z)\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})\\

把上两式加权和,有:

αf(y)+(1α)f(x)f(z)T(zz)+f(z)=f(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})

证毕。

定理:设SRnS\subset R^n 是非空开凸集,f:SR1f:S\rightarrow R^1 二阶可导,则ffSS 上的凸函数的充要条件是ff 的 Hessian 矩阵在SS 上都是半正定的。

当 Hessian 矩阵是正定的时,ff 是严格凸函数,但反之不成立。

若一个 MP 问题的目标函数ff 是凸函数,可行集XX 是凸集,则称规划问题为凸规划

定理:对于 MP 问题,若f,gif,g_i 函数均为凸函数,hjh_j 为线性函数,那么原问题是一个凸规划。

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

定理:凸规划的任一局部最优解都是它的整体最优解。

证明:设x\mathbf{x}^*MP 的一个局部最优解,故有一个邻域Nδ(x)N_\delta(\mathbf{x}^*) 满足:

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

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

f(αxˉ+(1α)x)αf(xˉ)+(1α)f(x)<αf(x)+(1α)f(x)=f(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 充分小时,αxˉ+(1α)x\alpha\bar{\mathbf{x}}+(1-\alpha)\mathbf{x}^* 可以成为邻域Nδ(x)N_\delta(\mathbf{x}^*) 中的点,故与x\mathbf{x}^* 是局部极小值相矛盾。证毕

# 一维搜索方法

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

方法 1:三分法。
给定一个搜索区间[a,b][a,b],我们可以在其中选择两个点x1,y1x_1,y_1,满足a<x1<y1<ba<x_1<y_1<b
5

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

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

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


方法 2:黄金分割法。

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

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

为保持这个过程可以永远进行下去,我们希望x1,y1x_1,y_1 在搜索区间[a,b][a,b] 中的位置和比例关系和y1,y2y_1,y_2 在搜索区间[x1,b][x_1,b]x,yx,y 的选择是对称的(因为也可能扔掉右边的区间):
7

故有比例关系:

1y1aba=y1x1bx11-\frac{y_1-a}{b-a}=\frac{y_1-x_1}{b-x_1}

设搜索区间[a,b][a,b] 长度为单位 1,设x1a=tx_1-a=t 那么by2=tb-y_2=ty1x1=12ty_1-x_1=1-2t。那么代入就会有:

11t1=12t1t1-\frac{1-t}{1}=\frac{1-2t}{1-t}

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

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

方法 3: Newton 法

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

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

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

x,f(x)f(tk)+f(tk)(xtk)+f(tk)2(xtk)2\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) 的极值点,就是让f(x)=0f'(x)=0,则有:

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

解得:

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

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

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

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

# 无约束最优化方法

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

定理:设f:RnR1f:R^n\rightarrow R^1xˉRn\bar{\mathbf{x}}\in R^n 处可微,若存在pRn\mathbf{p}\in R^n 使得:

f(xˉ)Tp<0\nabla f(\bar{\mathbf{x}})^T\mathbf{p}<0

则向量p\mathbf{p}ffxˉ\bar{\mathbf{x}} 处的下降方向。

证明:泰勒展开,得

f(xˉ+tp)=f(xˉ)+tf(xˉ)Tp+o(tp)f(\bar{\mathbf{x}}+t\mathbf{p})=f(\bar{\mathbf{x}})+t\nabla f(\bar{\mathbf{x}})^T\mathbf{p}+o(\|t\mathbf{p}\|)

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

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

故存在tt,使得f(xˉ+tp)<f(xˉ)f(\bar{\mathbf{x}}+t\mathbf{p})<f(\bar{\mathbf{x}})。证毕

定理:设f:RnR1f:R^n\rightarrow R^1 在点xRn\mathbf{x}^*\in R^n 处的 Hessian 矩阵2f(x)\nabla^2f(\mathbf{x}^*) 存在,若:

f(x)=0\nabla f(\mathbf{x}^*)=02f(x)\nabla^2f(\mathbf{x}^*) 正定

x\mathbf{x}^*UMP 的严格局部最优点。

证明:当 Hessian 矩阵正定时,说明ff 是一个可微严格凸函数,所以有:

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

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


最速下降法

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

f(xˉ)f(xˉ+tp)=tf(xˉ)Tpo(tp)f(\bar{\mathbf{x}})-f(\bar{\mathbf{x}}+t\mathbf{p})=-t\nabla f(\bar{\mathbf{x}})^T\mathbf{p}-o(\|t\mathbf{p}\|)

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

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

时,函数值下降最大。


共轭方向法

定义:设AAnn 阶实对称矩阵,对于非零向量p,qRn\mathbf{p},\mathbf{q}\in R^n,若有:

pTAq=0\mathbf{p}^TA\mathbf{q}=0

则称p,q\mathbf{p},\mathbf{q} 是相互AA 共轭的。对于非零向量组,满足:

piTApj=0,i,j=0,1,...,n1\mathbf{p}_i^TA\mathbf{p}_j=0,i,j=0,1,...,n-1

则称pi\mathbf{p}_iAA 的共轭方向组,称他们为一组AA 的共轭方向。

定理:设AAnn 阶实对称正定矩阵,pi\mathbf{p}_i 是非零向量,若p0,...,pn1\mathbf{p}_0,...,\mathbf{p}_{n-1} 是一组AA 共轭方向,则它们一定是线性无关的。

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

i=0n1αipi=0\sum_{i=0}^{n-1}\alpha_i \mathbf{p}_i=0

pjTA\mathbf{p}_j^TA 左乘上式:

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

由于AA 是正定的,且pi\mathbf{p}_i 是共轭方向,故当且仅当i=ji=j 时,有:

αipiTApi=0\alpha_i \mathbf{p}_i^TA\mathbf{p}_i=0

其他情况pjTApi=0\mathbf{p}_j^TA\mathbf{p}_i=0。又因为AA 正定,故piTApi>0\mathbf{p}_i^TA\mathbf{p}_i>0,故αi=0\alpha_i=0。同理可推出所有的αi=0\alpha_i=0。证毕。

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

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

其中AAnn 阶实对称正定矩阵。

定理:若p0,...,pn1\mathbf{p}_0,...,\mathbf{p}_{n-1} 为任意一组AA 的共轭方向,则由任意起始点x0\mathbf{x}_0 出发,依次沿着p0,...,pn1\mathbf{p}_0,...,\mathbf{p}_{n-1} 进行一维搜索,则最多经过nn 轮迭代可达整体最优解。

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

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

λ0=p0TAf(x1)p0TAp0\lambda_0=\frac{\mathbf{p}_0^TA\nabla f(\mathbf{x}_1)}{\mathbf{p}_0^TA\mathbf{p}_0}

  1. 继续在p1\mathbf{p}_1 方向上一维搜索。实际上,通式是:

pk+1=f(xk+1)+λkpkλk=pkTAf(xk+1)pkTApk\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}

为了简化记忆,实际上:

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

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

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

λk=f(xk+1)T[f(xk+1)f(xk)]f(xk)2\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={xRn    gi(x)0,hj(x)=0},i=1,...,p,j=1,...,qX=\{\mathbf{x}\in R^n\;|\;g_i(\mathbf{x})\leq 0,h_j(\mathbf{x})=0\},i=1,...,p,j=1,...,q

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

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

定理:设ffgi,iI(x)g_i,i\in I(\mathbf{x}^*)x\mathbf{x}^* 处可微,而gi,iI(x)g_i,i\notin I(\mathbf{x}^*)x\mathbf{x}^* 处连续,且hjh_jx\mathbf{x}^* 都连续可微。

并且,gi(x),iI(x),hj(x)\nabla g_i(\mathbf{x}^*),i\in I(\mathbf{x}^*),\nabla h_j(\mathbf{x}^*) 线性无关。若x\mathbf{x}^*MP 的局部最优解,则存在两组实数λi,μj\lambda_i^*,\mu_j^* 使得:

f(x)+iI(x)λigi(x)+jμjhj(x)=0λi0,iI(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\\ \lambda_i^*\geq 0,i\in I(\mathbf{x}^*)

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

式子f(x)+iI(x)λigi(x)+jμjhj(x)=0\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 条件的几何意义。首先若不考虑hjh_j 的部分,则有:

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

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

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

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

相比之下,若考虑hjh_j 的部分,则有:

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

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

若在 K-T 条件中,要求gi(x),iIg_i(\mathbf{x}),i\in Ix\mathbf{x}^* 处均可微,那么 K-T 条件可进一步为:

f(x)+i=1pλigi(x)+jμjhj(x)=0λi0,i=1,...,pλigi(x)=0,i=1,...,p\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(x)I(\mathbf{x}^*) 了呢?就是因为有了互补松紧条件 (λigi(x)=0\lambda_i^* g_i(\mathbf{x}^*)=0),若iI(xi\notin I(\mathbf{x}^*),则λi=0\lambda_i^*=0

引入 Lagrange 函数:

L(x,λ,μ)=f(x)+i=1pλigi(x)+j=1qμjhj(x)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 条件可写为:

{xL(x,λ,μ)=0λigi(x)=0,i=1,..,pλi0,i=1,...,p\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 乘子

定理:对于 MP 问题,若f,gi(iI),hjf,g_i(i\in I),h_jx\mathbf{x}^* 处连续可微,且x\mathbf{x}^* 满足 K-T 条件,且f,gi(iI),hjf,g_i(i\in I),h_jx\mathbf{x}^* 是凸函数,且hjh_j 是线性函数,则x\mathbf{x}^*MP 问题的整体最优解。

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

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

# 简约梯度法

考虑如下问题:

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

其中AA 是一个秩为mmm×nm\times n 的矩阵。现作如下约束非退化假设

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

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

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

仍有BxB+NxN=bB\mathbf{x}_B+N\mathbf{x}_N=\mathbf{b},故ff 可以表示为关于xN\mathbf{x}_N 的函数:

F(xN)=f(xB,xN)=f(B1bB1NxN,xN)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)

对它求梯度,有:

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

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


下面我们考虑如何构造搜索方向p=(pBpN)\mathbf{p}=\begin{pmatrix}\mathbf{p}_B\\ \mathbf{p}_N\end{pmatrix}, 对于一个x\mathbf{x} 的最大的mm 个分量组成的向量为xB>0\mathbf{x}_B>0(非退化假设),并记下标集合为IBI_B

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

我们这样选取pN\mathbf{p}_N 的分量:

iIBpi={F(xN)iF(xN)i0xiF(xN)iF(xN)i>0\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}

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

当然另一方面我们还希望保持Ax=bA\mathbf{x}=\mathbf{b} 的可行性,即:

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

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

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

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

{pN={F(xN)iF(xN)i0xiF(xN)iF(xN)i>0pB=B1NpN\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}

定理:对于非线性规划问题,设ff 是可微函数,p\mathbf{p} 是上述方法确定的搜索方向,则:

  1. p0\mathbf{p}\neq 0 时,p\mathbf{p}ffx\mathbf{x} 处的可行下降方向。
  2. p=0\mathbf{p}=0 得充要条件是x\mathbf{x} 是原问题的 K-T 点。

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

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

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

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

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

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

Bf(x)λB+BTμ=0\nabla_B f(\mathbf{x})-\lambda_B^*+B^T\mu^*=0

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

再根据简约梯度的式子,有:

f(x)λ+ATμ=0Bf(x)λB+BTμ=0λB=0μ=(BT)1Bf(x)f(x)λ+ATμ=0Nf(x)λN+NTμ=0λN=F(xN)\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}

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

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

证毕


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

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

于是tt 的上界为:

{+p0min1in{xipi    pi<0}otherwise\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. 选取初始可行点x\mathbf{x}
  2. IBI_Bx\mathbf{x} 的最大的mm 个分量的下标集合,A=(B,N)A=(B,N)。计算:

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

  1. 计算简约梯度:

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

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

{pN={F(xN)iF(xN)i0xiF(xN)iF(xN)i>0pB=B1NpN\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. 进行一维搜索,确定min  f(x+tp)min\;f(\mathbf{x}+t\mathbf{p})
  2. 重复 2-5 直到满足停止条件。

# 惩罚函数法

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

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

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

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

pc(c)=ci=1pmax(gi(x),0)2+c2j=1qhj(x)2p_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,gi,hif,g_i,h_i 均连续可微时,上述罚函数也是连续可微的。

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

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

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

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

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

Edited on Views times