第一章 基础知识

算法是一个满足下列条件的计算:

  • 输入:有满足给定约束条件的输入
  • 输出:能够输出满足给定约束条件的结果
  • 有穷性:有限步内必须停止
  • 确定性:每一步都是严格定义和确定的动作
  • 可行性:每一个动作都能够被精确地机械的执行

一般用算法的时间复杂度来度量算法的效率。

算法在机器上真正运行的时间取决于硬件性能,所以一般用该算法解决某问题所需的基本运算次数来表示算法的效率,这个次数通常还与问题的规模nn有关。所以时间复杂度一般都表示为输入规模的函数T(n)T(n)。

最后,相同规模的数据,也可能因为数据特点的不同,导致不同的基本运算次数。通常,算法的时间复杂度分为最坏情况下的时间复杂度W(n)W(n),以及平均情况下的时间复杂度A(n)A(n)。

A(n)=∑I∈SPItIA(n)=\sum_{I\in S}P_I t_I,其中SS是规模为nn的实例集,某个实例I∈SI\in S的概率为PIP_I,算法对实例II所需的基本运算次数是tIt_I。


设f,gf, g是定义域为自然数集N\N的函数

  1. 若存在正数c,n0c, n_0使得对于一切n⩾n0n\geqslant n_0,有0⩽f(n)⩽c⋅g(n)0\leqslant f(n)\leqslant c\cdot g(n)成立,则称f(n)f(n)的渐进的上界是g(n)g(n),记作f(n)=O(g(n))f(n)=O(g(n))
  2. 若存在正数c,n0c, n_0使得对于一切n⩾n0n\geqslant n_0,有0⩽c⋅g(n)⩽f(n)0\leqslant c\cdot g(n)\leqslant f(n)成立,则称f(n)f(n)的渐进的下界是g(n)g(n),记作f(n)=Ω(g(n))f(n)=\Omega(g(n))
  3. 若对于任意的正数cc都存在n0n_0,使得n⩾n0n\geqslant n_0时,有0⩽f(n)<c⋅g(n)0\leqslant f(n)< c\cdot g(n)成立,记作f(n)=o(g(n))f(n)=o(g(n))
  4. 若对于任意的正数cc都存在n0n_0,使得n⩾n0n\geqslant n_0时,有0⩽c⋅g(n)<f(n)0\leqslant c\cdot g(n)< f(n)成立,记作f(n)=ω(g(n))f(n)=\omega (g(n))
  5. 若f(n)=O(g(n))f(n)=O(g(n))且f(n)=Ω(g(n))f(n)=\Omega(g(n)),则g(n)g(n)称为渐进的紧的界,记为f(n)=Θ(g(n))f(n)=\Theta(g(n))

大OO记号的运算规则:

  • O(f)+O(g)=O(max⁡(f,g))O(f)+O(g)=O(\max(f, g))
  • O(f)+O(g)=O(f+g)O(f)+O(g)=O(f+g)
  • O(f)O(g)=O(fg)O(f)O(g)=O(fg)
  • O(cf(n))=O(f(n))O(cf(n))=O(f(n)),其中c>0c>0是一个常数
  • f=O(f)f=O(f)
  • 如果g(n)=O(f(n))g(n)=O(f(n)),则O(f)+O(g)=O(f)O(f)+O(g)=O(f)

设f,gf, g是定义域为自然数集N\N的函数

  • 如果lim⁡n→∞f(n)/g(n)\lim_{n\rightarrow \infin}f(n)/g(n)存在且等于某个常数c>0c>0,则f(n)=Θ(g(n))f(n)=\Theta(g(n))
  • 如果lim⁡n→∞f(n)/g(n)=0\lim_{n\rightarrow \infin}f(n)/g(n)=0,则f(n)=o(g(n))f(n)=o(g(n))
  • 如果lim⁡n→∞f(n)/g(n)=+∞\lim_{n\rightarrow \infin}f(n)/g(n)=+\infin,则f(n)=ω(g(n))f(n)=\omega(g(n))

设函数f,g,hf, g, h的定义域为N\mathbb{N},则:

  • 若f=O(g),g=O(h)f=O(g), g=O(h),有f=O(h)f=O(h)
  • 若f=Ω(g),g=Ω(h)f=\Omega(g), g=\Omega(h),有f=Ω(h)f=\Omega(h)
  • 若f=Θ(g),g=Θ(h)f=\Theta(g), g=\Theta(h),有f=Θ(h)f=\Theta(h)

多项式函数:f(n)=a0+a1n+⋯+adndf(n)=a_0+a_1n+\dots+a_dn^d,其中ad≠0a_d\neq 0。有f(n)=Θ(nd)f(n)=\Theta(n^d)

对数函数:对每个b>1,α>0b>1, \alpha>0,有log⁡bn=o(nα)\log_b n=o(n^\alpha),也即任何幂函数都比对数函数的阶要高。还有log⁡kn=Θ(log⁡ln)\log_k n=\Theta(\log_l n),无论底数如何,对数函数都是同阶的。

指数函数:对每个r>1r>1和每个d>0d>0,有nd=o(rn)n^d=o(r^n),也即任何指数函数都比多项式函数增长得快

阶乘函数:由Stirling 公式

n!=2πn(ne)n(1+Θ(1n))n!=\sqrt{2\pi n}\left(\frac n e\right)^n\left(1+\Theta\left(\frac 1 n\right)\right)

可得n!=o(nn),n!=ω(2n),log⁡(n!)=Θ(nlog⁡n)n!=o(n^n), n!=\omega(2^n), \log(n!)=\Theta(n\log n)


调和级数,可以用积分做其渐进的界:

ln⁡(n+1)=1+∫1n1xdx⩾1+∑k=2n1k=∑k=1n1k⩾∫1n+11xdx=ln⁡n\ln (n+1)=1+\int_1^n \frac 1 x dx\geqslant 1+\sum_{k=2}^n \frac 1 k=\sum_{k=1}^n\frac 1 k \geqslant \int_{1}^{n+1} \frac 1 x dx=\ln n

所以:

∑k=1n1k=Θ(ln⁡n)\sum_{k=1}^n\frac 1 k=\Theta(\ln n)


主定理:设a⩾1,b>1a\geqslant 1, b>1为常数,f(n)f(n)为函数,T(n)T(n)为非负整数,且:

T(n)=aT(nb)+f(n)T(n)=aT(\frac n b)+f(n)

则:

  • 若f(n)=O(nlog⁡ba−ε),ε>0f(n)=O(n^{\log _ba-\varepsilon}), \varepsilon>0,那么T(n)=Θ(nlog⁡ba)T(n)=\Theta(n^{\log_ba})
  • 若f(n)=Θ(nlog⁡ba)f(n)=\Theta(n^{\log_ba}),那么T(n)=Θ(nlog⁡balog⁡n)T(n)=\Theta(n^{\log_b a}\log n)
  • 若f(n)=O(nlog⁡ba+ε),ε>0f(n)=O(n^{\log _ba+\varepsilon}), \varepsilon>0,且对于某个常数c<1c<1和所有充分大的nn都有af(n/b)⩽cf(n)af(n/b)\leqslant cf(n),那么T(n)=Θ(f(n))T(n)=\Theta(f(n))

第二章 分治算法

1、基本概念

分治算法基本思想是将一个规模为nn的问题以某种方式分解为kk个规模较小的子问题,这些子问题非常小以至于能在常数时间内解决,最后将这些子问题的解合并为原问题的解。这三步就是分治算法的三个步骤:

  1. 分(Divide):将大规模问题分割成若干个更小规模的子问题。如果子问题的规模不够小,则再继续划分,如此递归地进行下去。
  2. 治(Conquer):解决这些子问题。
  3. 合(Combine):将这些子问题的解合并为原问题的解。

其中分这个步骤是分治算法基础和关键,一般要遵循两个原则:

  • 平衡子问题原则:分割出的若干个子问题,其规模最好大致相当
  • 独立子问题原则:分割出的若干个子问题,之间的重叠越少越好,最好是互相独立的

一般地,分治算法时间复杂度可以写为以下递推形式:

W(n)=W(∣P1∣)+W(∣P2∣)+⋯+W(∣Pk∣)+f(n)W(c)=C\begin{align*} W(n)=W(|P_1|)+W(|P_2|)+\dots+W(|P_k|)+f(n) \\ W(c)=C \end{align*}

其中∣Pi∣|P_i|是第ii个子问题的规模,f(n)f(n)是合并子问题解的时间开销,CC是直接求解规模为cc的子问题的时间开销。


根据分解出的子问题规模,具体还有以下两种常见的递推形式:

T(n)=∑i=1kaiT(n−i)+f(n)T(n)=aT(n/b)+d(n)\begin{align*} T(n)=\sum_{i=1}^k a_i T(n-i)+f(n) \\ T(n)=aT(n/b)+d(n) \end{align*}

例如:

  • 汉诺塔问题T(n)=2T(n−1)+1T(n)=2T(n-1)+1
  • 二分查找问题T(n)=T(n/2)+1T(n)=T(n/2)+1
  • 归并排序问题T(n)=2T(n/2)+nT(n)=2T(n/2)+n。

对于第一个方程,可以用迭代法、递归树、尝试法求解;对于第二个方程,易看出这是主定理的形式。

当d(n)d(n)为常数时,由主定理有:

T(n)={Θ(nlog⁡ba)a≠1Θ(log⁡n)a=1T(n)=\begin{cases} \Theta(n^{\log_ba}) & a\neq 1\\ \Theta(\log n) & a=1 \end{cases}

当d(n)=cnd(n)=cn时,分别对应主定理的三种情况:

T(n)={Θ(n)a<bΘ(nlog⁡n)a=bΘ(nlog⁡ba)a>bT(n)=\begin{cases} \Theta(n) & a<b\\ \Theta(n\log n) & a=b\\ \Theta(n^{\log_ba}) & a>b \end{cases}

2、实例

2.1、逆序对问题

给定一个包含nn个元素的数组AA,一个逆序对是一个满足i<ji<j且A[i]>A[j]A[i]>A[j]的有序对。求逆序对的数量。

  • 分:当数组规模为n>2n>2时,将数组分为两个规模为n/2n/2的子数组
  • 治:递归地求解两个子数组的逆序对数量。如果n=1n=1,则逆序对数量为 0;如果n=2n=2,则逆序对数量为 0 或 1。
  • 合:原数组的逆序对数量=两个子数组之间的逆序对数量+跨两个子数组的逆序对数量。

治的时候同时使元素按升序排列,并在并的时候按照升序合并,这样可以使计算跨两个子数组的逆序对的复杂度降为O(n)O(n)。(可见,这就是归并排序。所以求逆序对可以再归并排序的基础上进行)


时间复杂度递推式:

W(n)=2W(n/2)+O(n)W(1)=0,W(2)=1\begin{align*} W(n)=2W(n/2)+O(n)\\ W(1)=0, W(2)=1 \end{align*}

所以W(n)=O(nlog⁡n)W(n)=O(n\log n)

2.2、芯片测试问题

现有nn个芯片,其中好芯片至少比坏芯片多一片。每次拿两个芯片测试,每个芯片会报告另一个芯片是好或者坏。好芯片的报告总是正确的,坏芯片的报告可能是正确的,也可能是错误的。设计一个算法,使用最少的测试次数找到一个好芯片。

考虑用其他n−1n-1个芯片对剩下的芯片 A 进行测试。首先,好芯片的数量至少为⌊n2⌋+1\lfloor\frac{n}{2}\rfloor+1,这些好芯片对 A 的报告一定是正确的

  • A 是好的,则至少有⌊n2⌋\lfloor\frac{n}{2}\rfloor个芯片报告 A 为好
  • A 是坏的,则至少有⌊n2⌋+1\lfloor\frac{n}{2}\rfloor+1个芯片报告 A 为坏

另外,两个芯片都报告对方是好时,两个芯片都好或者都坏。


于是可以考虑如下分治策略:

  • 分:当n>3n>3时,两两一组做测试,互相报告为好的组,任留一片,另一片丢弃;其他情况,两片都丢弃。当nn为奇数时,最后剩下的一个芯片用其余n−1n-1个芯片测试,由上面的结论,可以直接判断出该芯片是好还是坏。如果是好的,算法结束;如果是坏的,直接丢弃。
  • 治:n<=3n<=3时,一次测试即可确定好芯片。
  • 合:治的时候,已经得到了

下面证明算法正确性。假设nn为偶数,考虑分的时候,A、B 都好的有ii组,一好一坏的有jj组,A、B 都坏的有kk组。则有:

2i+2j+2k=n2i+j>2k+j\begin{align*} 2i+2j+2k=n\\ 2i+j>2k+j \end{align*}

淘汰后,好芯片数量为ii,坏芯片数量为kk,由上面的式子,有i>ki>k。所以“好芯片至少比坏芯片多一片”的性质始终保留。nn为奇数时,通过轮空处理,如果算法没结束,丢弃的一定是坏芯片,性质仍然保留,转化为nn为偶数的情况。

最终剩余 3 片以内时:

  • n=1,2n=1, 2,剩下的芯片一定都是好的
  • n=3n=3,好芯片至少两片。任取两个芯片,做一次测试,一定可以判断出三个芯片的好坏情况。

时间复杂度递推式:

W(n)=W(n/2)+O(n)W(3)=1,W(2)=W(1)=0\begin{align*} W(n)=W(n/2)+O(n)\\ W(3)=1, W(2)=W(1)=0 \end{align*}

最终有W(n)=O(n)W(n)=O(n)

2.3、快速排序

给定一个长度为nn的数组 A,输出排序后的数组。

  • 分:当nn>1 时,选取某个元素xx为基准,将数组分为两个子数组 A1、A2,使得 A1 中的元素都小于等于xx,A2 中的元素都大于xx。接下来递归地对 A1、A2 进行排序。
  • 治:当n<=1n<=1时,数组已排好序。
  • 合:合并 A1、x、A2。由于子数组都是排好序的,[A1, x, A2]就是有序的。

快速排序的时间复杂度与选择的基准元素有关,划分出的两个子数组的规模直接影响了排序效率。

最坏的情况下,每次一个子数组长度满的,另一个子数组空的:

W(n)=W(n−1)+n−1W(1)=W(0)=0\begin{align*} W(n)=W(n-1)+n-1\\ W(1)=W(0)=0 \end{align*}

所以W(n)=n(n−1)/2W(n)=n(n-1)/2

最好的情况下,每次划分都能均匀划分:

T(n)=2T(n/2)+n−1T(1)=T(0)=0\begin{align*} T(n)=2T(n/2)+n-1\\ T(1)=T(0)=0 \end{align*}

所以T(n)=Θ(nlog⁡n)T(n)=\Theta(n\log n)


考虑计算平均时间复杂度。假设每次选择首元素划分后,首元素位于第ii个位置(从 1 开始),那么有:

T(n)=T(i−1)+T(n−i)+n−1T(n)=T(i-1)+T(n-i)+n-1

首元素处于每个位置的概率都是1/n1/n,所以:

A(n)=1n∑i=1n[T(i−1)+T(n−i)+n−1]A(1)=A(0)=0\begin{align*} A(n)=\frac 1 n\sum_{i=1}^n[T(i-1)+T(n-i)+n-1]\\ A(1)=A(0)=0 \end{align*}

可算出T(n)=Θ(nlog⁡n)T(n)=\Theta(n\log n)

2.4、快速幂

给定整数a,ba, b,求aba^b。

  • 分:b>=2b>=2时,对指数分解。当bb为偶数时,考虑计算ab/2a^{b/2};当bb为奇数时,考虑计算a(b−1)/2a^{(b-1)/2}
  • 治:b=1b=1,aa就是子问题结果
  • 和:对子问题结果平方,再视bb的奇偶性决定是否再乘以aa

时间复杂度的递推式:

W(n)=W(n/2)+Θ(1)W(1)=0\begin{align*} W(n)=W(n/2)+\Theta(1)\\ W(1)=0 \end{align*}

可得W(n)=Θ(log⁡n)W(n)=\Theta(\log n)


利用相同的思想可以得到快速矩阵幂。例如斐波那契数列的计算可以写成矩阵形式:

[Fn+2Fn+1Fn+1Fn]=[Fn+1FnFnFn−1][1110]=[1110]n[1110]=[1110]n+1\begin{bmatrix} F_{n+2} & F_{n+1}\\ F_{n+1} & F_n \end{bmatrix}=\begin{bmatrix} F_{n+1} & F_n\\ F_n & F_{n-1} \end{bmatrix}\begin{bmatrix} 1 & 1\\ 1 & 0 \end{bmatrix}=\begin{bmatrix} 1 & 1\\ 1 & 0 \end{bmatrix}^n\begin{bmatrix} 1 & 1\\ 1 & 0 \end{bmatrix}=\begin{bmatrix} 1 & 1\\ 1 & 0 \end{bmatrix}^{n+1}

接下来可以用矩阵快速幂的方法计算。

2.5、选择问题

给定一个长度为nn的数组 A ,找最大和最小的元素。

  • 直接遍历。最坏的时间复杂度为W(n)=n−1+n−2=n−3W(n)=n-1+n-2=n-3

  • 分治算法。两两一组,可以分为⌊n/2⌋\lfloor n/2\rfloor组。组内比较大小,得到⌈n/2⌉\lceil n/2\rceil个“较小”/“较大”(奇数时为⌊n/2⌋+1\lfloor n/2\rfloor+1)。然后遍历这些“较小”/“较大”,得到最小/最大的元素,还需要比较2(⌈n/2⌉−1)2(\lceil n/2\rceil-1)次。所以最坏时间复杂度为W(n)=⌊n/2⌋+2⌈n/2⌉−2=⌈3n/2⌉−2W(n)=\lfloor n/2\rfloor+2\lceil n/2\rceil-2=\lceil 3n/2\rceil-2,比直接遍历要好。可以证明,这就是时间复杂度最低的算法。


给定一个长度为nn的数组 A ,找第二大的元素。

  • 直接遍历两次。最坏的时间复杂度为W(n)=n−1+n−2=2n−3W(n)=n-1+n-2=2n-3

  • 锦标赛算法。每次两两一组比赛,胜者进入下一轮。由于第二大的元素只可能被最大的元素淘汰。所以可以每次做比较时,将被败者元素添加到胜者元素对应的链表上。最后查找最大元素的链表中最大的元素即可。

    最大元素能一直进入下一轮,比较⌈log⁡2n⌉\lceil \log_2 n\rceil次,在链表中查找最大元素需要比较n−1n-1次,所以最终其链表长度为⌈log⁡2n⌉−1\lceil \log_2 n\rceil-1。每轮比赛每组都淘汰一个元素,最终一共淘汰n−1n-1个元素,进行了n−1n-1次比较。所以时间复杂度为W(n)=⌈log⁡2n⌉−1+n−1=n+⌈log⁡2n⌉−2W(n)=\lceil \log_2 n\rceil-1+n-1=n+\lceil \log_2 n\rceil-2。可以证明,这就是时间复杂度最低的算法。


给定一个长度为nn的数组 A ,找第 k 大的元素。这是最一般的选择问题。

  • 直接排序后查找。时间复杂度为O(nlog⁡n)O(n\log n)

  • 分治算法。每次选取一个基准元素m∗m^*,比其小的划分为一个子数组,比其大的划分为另一个子数组。根据子数组的元素个数,决定进入哪一个子数组进行递归。

    如果选择算法的时间复杂度为T(n)T(n),则选择基准元素的复杂度应该为T(cn)T(cn),其中c<1c<1。最坏情况下,每次递归调用都进入规模较大的子数组,解决子问题的时间复杂度应该为T(dn)T(dn),其中d>1d>1。而且要有c+d<1c+d<1,才能使最终复杂度达到O(n)O(n)。

    考虑五个元素一组,找到每组的中位数,然后找到所有中位数的中位数,作为基准元素。然后将整个数组按下图分为四块:

    select-algorithm

    每列就是刚才的一组,设其是从大到小排序的,所以组的中位数恰好都在第三行。列之间根据中位数大小排序基准元素就在第三行的中间。显然 C 中元素都比m∗m^*小,B 中元素都比m∗m^*大。需要遍历 A、D 两块来确定其中哪些元素大于/小于m∗m^*。这样,得到划分后的子数组。然后递归调用。

    不妨设nn是 5 的倍数,且n/5=2r+1n/5=2r+1为奇数。有∣A∣=∣D∣=2r,∣B∣=∣C∣=3r+2|A|=|D|=2r, |B|=|C|=3r+2。最坏的情况下,递归调用进入的子数组规模为∣A∣+∣D∣+∣C∣=7r+2<0.7n|A|+|D|+|C|=7r+2<0.7n。

    所以最坏情况下的时间复杂度W(n)⩽W(n/5)+W(7n/10)+tnW(n)\leqslant W(n/5)+W(7n/10)+tn。由主定理,有W(n)=O(n)W(n)=O(n)。不等号右侧第一项对应调用选择算法查找中位数的中位数,第二项对应递归调用,第三项对应找每组中位数以及处理 A、D 两块的时间。

2.6、多项式在 1 的全体2n2n次方根的值

1 在复数域上开2n2n次方,有2n2n个根ωi=cos⁡πjn+isin⁡πjn\omega_i=\cos \frac{\pi j}{n}+i\sin \frac{\pi j}{n},i=0,1,…,2n−1i=0, 1, \dots, 2n-1。给定一个多项式A(x)=a0+a1x+⋯+an−1xn−1A(x)=a_0+a_1x+\dots+a_{n-1}x^{n-1},求所有A(ωi)A(\omega_i)。

  • 根据定义,一个个的求。求单个A(ωi)A(\omega_i)的时间复杂度为O(n2)O(n^2),总的时间复杂度为O(n3)O(n^3)

  • 迭代法。有迭代式Ai(x)=an−i+xAi−1(x),A1(x)=an−1A_i(x)=a_{n-i}+xA_{i-1}(x), A_1(x)=a_{n-1},所以An(x)=A(x)A_n(x)=A(x)。
    求单个值的时间复杂度为O(n)O(n),总的时间复杂度为O(n2)O(n^2)

  • 分治算法。不妨设n=2rn=2r,考虑多项式A0(x)=a0+a2x+⋯+a2rx,A1(x)=a1+a3x+⋯+a2r−1(x)A_0(x)=a_0+a_2 x+\cdots+a_{2r}x, A_1(x)=a_1+a_3 x+\cdots+a_{2r-1}(x)。则A(x)=A0(x2)+xA1(x2)A(x)=A_0(x^2)+xA_1(x^2)。所以可以递归地求解。

    有ωi2=ω(2i)%(2n)\omega_i^2=\omega_{(2i)\%(2n)}。所以A(ωi)=A0(ω(2i)%(2n))+ωiA1(ω(2i)%(2n)),∀jA(\omega_i)=A_0(\omega_{(2i)\%(2n)})+\omega_i A_1(\omega_{(2i)\%(2n)}), \forall j,也即原问题可以划分为两个规模减半的子问题。所以整个问题时间复杂度为T(n)=2T(n/2)+O(n)T(n)=2T(n/2)+O(n),由主定理有T(n)=O(nlog⁡n)T(n)=O(n\log n)。

2.7、平面点集凸包问题

给定一个平面上的点集,找到包含所有点的最小凸多边形。

考虑分治算法。首先找到纵坐标最大和最小的两个点,用它们之间的连线dd将点集划分为左右两部分。然后找到距离dd最远的点PP,这个点与dd的两个端点连线分别为a,ba, b,则a,b,da, b, d构成一个三角形。三角形内的点直接删去,a 及其外侧的点构成一个新的点集,b 及其外侧的点构成另一个新的点集。递归地求解。

初始划分的时间复杂度为O(n)O(n)。每次根据dd找到PP的时间复杂度为O(n)O(n),然后划分子问题的时间复杂度为O(n)O(n)。最坏的情况下,每次划分出的子问题规模为n−1n-1,时间复杂度为W(n)=W(n−1)+O(n)W(n)=W(n-1)+O(n),可得W(n)=O(n2)W(n)=O(n^2)。所以总的时间复杂度为O(n2)O(n^2)。

3、分治算法的改进

3.1、减少子问题个数

考虑分治算法时间复杂度递推式:

W(n)=aW(n/b)+d(n)W(n)=aW(n/b)+d(n)

当a>ba>b,d(n)d(n)不大时,由主定理,有W(n)=Θ(nlog⁡ba)W(n)=\Theta(n^{\log_ba}),此时减少aa可以降低W(n)W(n)的阶。

利用子问题之间的依赖关系,使得某些子问题的解可以通过组合其他子问题的解得到。这样,可以减少子问题的个数,降低时间复杂度。


考虑两个nn位的二进制数X,YX, Y相乘。直接相乘,需要O(n2)O(n^2)次乘法运算。

将每个数分为两部分X=X1⋅2n/2+X0,Y=Y1⋅2n/2+Y0X=X_1\cdot 2^{n/2}+X_0, Y=Y_1\cdot 2^{n/2}+Y_0,则其乘积可以写为:

X⋅Y=(X1⋅2n/2+X0)(Y1⋅2n/2+Y0)=X1Y1⋅2n+(X1Y0+X0Y1)⋅2n/2+X0Y0X\cdot Y=(X_1\cdot 2^{n/2}+X_0)(Y_1\cdot 2^{n/2}+Y_0)=X_1Y_1\cdot 2^n+(X_1Y_0+X_0Y_1)\cdot 2^{n/2}+X_0Y_0

时间复杂度W(n)=4W(n/2)+O(n)W(n)=4W(n/2)+O(n),由主定理,有W(n)=O(nlog⁡24)=O(n2)W(n)=O(n^{\log_2 4})=O(n^2)。

寻找子问题之间的依赖关系,可以发现:

X1Y0+X0Y1=(X1+X0)(Y1+Y0)−X1Y1−X0Y0X_1Y_0+X_0Y_1=(X_1+X_0)(Y_1+Y_0)-X_1Y_1-X_0Y_0

所以仅需要三次乘法运算,就可以得到X⋅YX\cdot Y。时间复杂度降为W(n)=3W(n/2)+O(n)W(n)=3W(n/2)+O(n),由主定理,有W(n)=O(nlog⁡23)=O(n1.59)W(n)=O(n^{\log_2 3})=O(n^{1.59})。


Strassen 矩阵乘法,也用了类似的思想。将两个矩阵相乘,考虑将每个矩阵分为四个子矩阵,可以得到:

[A11A12A21A22][B11B12B21B22]=[C11C12C21C22]\begin{bmatrix} A_{11} & A_{12}\\ A_{21} & A_{22} \end{bmatrix}\begin{bmatrix} B_{11} & B_{12}\\ B_{21} & B_{22} \end{bmatrix}= \begin{bmatrix} C_{11} & C_{12}\\ C_{21} & C_{22} \end{bmatrix}

如果只是普通的分块计算:

C11=A11B11+A12B21C12=A11B12+A12B22C21=A21B11+A22B21C22=A21B12+A22B22\begin{align*} C_{11}=A_{11}B_{11}+A_{12}B_{21}\\ C_{12}=A_{11}B_{12}+A_{12}B_{22}\\ C_{21}=A_{21}B_{11}+A_{22}B_{21}\\ C_{22}=A_{21}B_{12}+A_{22}B_{22} \end{align*}

时间复杂度为W(n)=8W(n/2)+O(n2)W(n)=8W(n/2)+O(n^2),由主定理,仍然是W(n)=O(n3)W(n)=O(n^3)。

考虑这样计算:

M1=A11(B12−B22)M2=(A11+A12)B22M3=(A21+A22)B11M4=A22(B21−B11)M5=(A11+A22)(B11+B22)M6=(A12−A22)(B21+B22)M7=(A11−A21)(B11+B12)\begin{align*} M_1=A_{11}(B_{12}-B_{22})\\ M_2=(A_{11}+A_{12})B_{22}\\ M_3=(A_{21}+A_{22})B_{11}\\ M_4=A_{22}(B_{21}-B_{11})\\ M_5=(A_{11}+A_{22})(B_{11}+B_{22})\\ M_6=(A_{12}-A_{22})(B_{21}+B_{22})\\ M_7=(A_{11}-A_{21})(B_{11}+B_{12}) \end{align*}

则:

C11=M5+M4−M2+M6C12=M1+M2C21=M3+M4C22=M5+M1−M3−M7\begin{align*} C_{11}=M_5+M_4-M_2+M_6\\ C_{12}=M_1+M_2\\ C_{21}=M_3+M_4\\ C_{22}=M_5+M_1-M_3-M_7 \end{align*}

所以时间复杂度为W(n)=7W(n/2)+O(n2)W(n)=7W(n/2)+O(n^2),由主定理,有W(n)=O(nlog⁡27)=O(n2.81)W(n)=O(n^{\log_2 7})=O(n^{2.81})。

目前最好的算法是 Coppersmith-Winograd 算法,时间复杂度为O(n2.376)O(n^{2.376})。

3.2、增加预处理

用平面最邻近点对问题说明:给定一个平面点集,找到之间距离最短的一对点。

朴素方法,两两枚举。共有Cn2C_n^2对点,时间复杂度为O(n2)O(n^2)。

分治方法,将点集分为两部分,分别求解,然后考虑跨两部分的最邻近点对,一共三部分。

  • 分割考虑做中垂线(也即按照 x 坐标分割),分为大小相近的两部分。一开始按 x 坐标进行排序O(nlog⁡n)O(n\log n),后续的分割只需要O(n)O(n)时间,总共是O(nlog⁡n)O(n\log n)
  • 对于找跨两部分的点对,设两部分中最邻近点对的距离分别是δ0,δ1\delta_0, \delta_1,则只需考虑中垂线两边δ=min⁡(δ0,δ1)\delta=\min(\delta_0, \delta_1)之内的点。假设中垂线左侧,距中垂线距离小于δ\delta的一个点(x1,y1)(x_1, y_1),只需要考虑中垂线右侧,同样到中垂线距离小于δ\delta,且纵坐标范围是[y1−δ,y1+δ][y_1-\delta, y_1+\delta]的点。可以证明,这样的点至多有 6 个。找到点(x1,y1)(x_1, y_1)需要O(n)O(n)时间,在 y 坐标有序的情况下,找到另一侧满足要求的点需要O(n)O(n)时间,检查是O(1)O(1)的,总共也是O(nlog⁡n)O(n\log n)的(按 y 进行排序)。

所以总的时间复杂度递推式为W(n)=2W(n/2)+O(nlog⁡n)W(n)=2W(n/2)+O(n\log n),用递归树可解得W(n)=O(nlog⁡2n)W(n)=O(n\log^2 n)


上面的分治方法中,每次递归都要调用一次排序(对 y 进行排序),这是因为分割的时候是对 x 进行排序,不能保证顶点数组中 y 的坐标也是有序的。有没有可能,分割的时候,同时能得到 x 的有序性和 y 的有序性,且还能够正常将每个点的 x 和 y 对应起来?

考虑预处理,一开始按 y 也进行排序,不过要带上额外信息,用于确定其属于那个点。

用 C++数据结构来说明:一开始有顶点坐标数组vector<pair<float, float>> points。创建pair<float, int> x,其中第一个元素是某个顶点的 x 坐标,第二个元素是顶点 ID(可以取在原顶点数组中的下标)。然后按 x 坐标排序,得到有序的vector<pair<float, int>>。再创建pair<float, int>,其中第一个元素是 y 坐标,第二个元素是顶点 ID 。按 y 坐标排序,得到有序的vector<pair<float, int>>。

然后进行分割。可以在O(n)O(n)的时间内得到vector<pair<float, int>> x_1, x_2。根据顶点 ID,扫描vector<pair<float, int>> y,顶点 ID 相同,说明是同一个点的 y 坐标,划分到对应的部分。可以在O(n)O(n)内得到vector<pair<float, int>> y_1, y_2,且仍然是有序的。

这样时间复杂度递推式为W(n)=2W(n/2)+O(n)W(n)=2W(n/2)+O(n),由主定理,有W(n)=O(nlog⁡n)W(n)=O(n\log n)。

第三章 动态规划

1、基本概念

多阶段决策问题是指,求解的问题求解的问题可以划分为一系列相互联系的阶段,在每个阶段都需要作出决策,且一个阶段决策的选择会影响下一个阶段的决策,从而影响整个过程的活动路线,求解的目标是选择各个阶段的决策使整个过程达到最优。

动态规划(Dynamic Programming)就是一种在研究多阶段决策问题时提出的一种解方法,其基本思想是把多阶段过程转化为一系列单阶段问题,然后逐个求解。动态规划常常用于求解以时间划分阶段的动态过程的优化问题,对于一些与时间无关的静态规划,也可以以人为引入时间因素,然后同样适用动态规划地方法求解。


  • 阶段:利用动态规划求解问题,需要将问题恰当地划分为若干个相互联系的阶段
  • 状态:每个阶段开始时,问题或者系统所处的客观状况。状态既是某个阶段的某个起点,也是前一个阶段的某个终点,一个阶段可以有若干种可能的状态。
  • 策略:每个阶段都需要作出决策,决策使得系统从一个状态转移到另一个状态。各个阶段的决策构成一个决策序列,这个序列就称之为一个策略。从某个阶段开始到终止阶段的过程称为一个子过程,对应的策略称为一个子策略。

状态的无后效性是指,某个阶段的状态给定之后,则该阶段之后的过程的发展不受该阶段以前各段状态的影响,也就是说状态具有马尔可夫性。适用于动态规划求解的问题,其中各个状态必须具有无后效性。


动态规划的实质是分治+减少冗余计算。

动态规划也是将原问题分解为若干个子问题,先求解子问题,然后根据子问题的解得到原问题的解。

与分治不同的是,动态规划中子问题往往不是相互独立的,会出现相同的子问题。如果使用分治方法求解,会有很多重复计算。动态规划用一个表来保存子问题的解,自底向上计算,最终求出原问题的解。


利用动态规划求解问题的一般步骤:

  1. 找出最优解的性质,并刻画其结构特征
  2. 递归地定义最优值,也即写出动态规划方程(状态转移方程)
  3. 自底向上计算最优值
  4. 根据计算最优值时得到的信息,构造最优解(可选)

2、Bellman 最优性原理

如果某个问题的最优策略的子策略总是最优的,则称该问题满足 Bellman 最优性原理。对于满足 Bellman 最优性原理的问题,如果其某个策略的子策略不是最优的,则该策略也不是最优的。


有向带权图GG中,从顶点viv_i到vjv_j的最短路径问题是满足最优性原理的。

证:假设该问题不满足最优性原理,则存在一条viv_i到vjv_j的最短路径vi→u→w→vjv_i\rightarrow u \rightarrow w \rightarrow v_j,其中的子路径u→w→vju \rightarrow w\rightarrow v_j不是uu到vjv_j的最短路径。

假设uu到vjv_j的最短路径是u→w′→vju\rightarrow w'\rightarrow v_j,则路径vi→u→w′→vjv_i\rightarrow u \rightarrow w' \rightarrow v_j比原来的路径更短,与原来的路径是最短路径的假设矛盾。

由反证法可知,从顶点viv_i到vjv_j的最短路径问题满足最优性原理。


无向无权图GG中,从顶点qq到tt的最长路径问题不满足最优性原理。设GG是一个环q↔r↔s↔t↔qq\leftrightarrow r \leftrightarrow s \leftrightarrow t \leftrightarrow q。

qq到tt的最长路径是q→r→tq\rightarrow r\rightarrow t。但是qq到rr的最长路径是q→s→t→rq\rightarrow s\rightarrow t\rightarrow r;rr到tt的最长路径是r→q→s→tr\rightarrow q\rightarrow s\rightarrow t。两个子问题的最优策略组合起来,不是整个问题的最优策略。说明该问题不满足最优性原理。


动态规划方法对问题的有效性,取决于问题本身是否满足:

  • 最优子结构:也称为优化原则,是指一个最优决策序列的任何子序列本身一定是相对于子序列初始和结束状态最优的决策序列。
  • 重叠子问题:递归求解时,会需要反复求解相同的子问题。动态规划方法会将子问题的解保存在一个表中,能够避免重复计算。

动态规划算法设计的要点:

  1. 问题要求什么?约束条件是什么?
  2. 如何划分子问题?
  3. 原问题的最优值与子问题的最优值之间的关系是什么?(状态转移方程)
  4. 是否满足最优子结构?
  5. 最小的子问题是什么?其解如何计算?(边界条件)

3、实例

3.1、多起点多终点的最短路径问题

给定一个有向带权图GG,起点集{S1,S2,…,Sn}\{S_1, S_2, \dots, S_n\},终点集{T1,T2,…,Tm}\{T_1, T_2, \dots, T_m\},求出起点在起点集,终点在终点集的最短路径。

shortest-path


蛮力算法,穷举每一条可能的路径。假定起点到终点的(平均)要经过kk条边,则时间复杂度达到O(n2k)O(n2^k)。

动态规划算法,考虑从终点开始,逐步前推。先求出起点集为{Ci}\{C_i\},终点为{Tj}\{T_j\}的最短路径F(Ci)=min⁡j{CiTj}F(C_i)=\min_{j}\{C_iT_j\}。然后求出起点集为{Bk}\{B_k\},终点为{Tj}\{T_j\}的最短路径F(Bk)=min⁡j{BkCj+F(Cj)}F(B_k)=\min_{j}\{B_kC_j+F(C_j)\}。依次类推,最终有F(Sl)=min⁡m{SlAm+F(Am)}F(S_l)=\min_{m}\{S_lA_m+F(A_m)\}。最小的F(Sl)F(S_l)对应的路径即为所求。


能使用动态规划算法,是因为满足最优子结构。即全局最短路径的子路径,也一定是相对这个子路径起点和终点的最短路径。

3.2、矩阵链相乘

给定矩阵序列A1,A2,…,An\boldsymbol{A}_1, \boldsymbol{A}_2, \dots, \boldsymbol{A}_n,其中Ai\boldsymbol{A}_i的规模为Pi−1×PiP_{i-1}\times P_i,求出最优的矩阵相乘顺序,使得计算元素相乘的总次数最少。(矩阵的行数和列数限定整个序列的顺序不变,利用结合律,通过加括号的方法,得到不同的计算次数)


蛮力算法,考虑穷举每一种可能的加括号方式。加完括号的序列可以写成一个二叉树,树的每个叶子节点都对应一个矩阵,每个子树对应着一个一对括号。由叶子节点开始,逐步向上计算,最终得到根节点对应的矩阵。根节点矩阵就对应着答案,整个过程就对应着计算的过程。

设nn个叶子节点的二叉树有xnx_n种,有递推公式:

xn+1=∑i=1nxixn+1−i,x1=1,x2=1x_{n+1}=\sum_{i=1}^n x_i x_{n+1-i}, x_1=1, x_2=1

也即ii个叶子结点的二叉树和n+1−in+1-i个叶子节点的二叉树,分别作为根节点的两个子树,得到n+1n+1个叶子节点的二叉树。这个形式与卡特兰数(Catalan Number)的定义本质上相同:

Cn=∑i=0n−1CiCn−1−i,C0=1,C1=1C_{n}=\sum_{i=0}^{n-1}C_{i}C_{n-1-i}, C_0=1, C_1=1

也即有xn+1=Cnx_{n+1}=C_n。

所以对于长度为n+1n+1的矩阵序列,其可能的运算顺序有CnC_n种,利用 Stirling 公式,可得穷举的时间复杂度为:

W(n)=Ω(Cn)=Ω(1n+1(2n)!n!n!)=Ω(4nn3/2)W(n)=\Omega(C_n)=\Omega(\frac{1}{n+1}\frac{(2n)!}{n!n!})=\Omega\left(\frac{4^n}{n^{3/2}}\right)

这是一个指数级别的时间复杂度。


动态规划算法。考虑某个矩阵链AiAi+1…Aj\boldsymbol{A}_i\boldsymbol{A}_{i+1}\dots\boldsymbol{A}_j,记其最少运算次数为m[i,j]m[i, j]。假设其最后一次相乘是在kk处(运算树的根节点对应的位置),也即最后相当于是AiAi+1…Ak\boldsymbol{A}_i\boldsymbol{A}_{i+1}\dots\boldsymbol{A}_k和Ak+1Ak+2…Aj\boldsymbol{A}_{k+1}\boldsymbol{A}_{k+2}\dots\boldsymbol{A}_j相乘。假设两个子部分都采用了最优的运算次序,则这种情况下的乘法次数应该是:

m[i,k]+m[k+1,j]+Pi−1PkPjm[i, k]+m[k+1, j]+P_{i-1}P_kP_j

所以每一个问题,可以遍历最后一次乘法出现的位置,写出递推式如下:

m[i,j]={0,i=jmin⁡i⩽k<j{m[i,k]+m[k+1,j]+Pi−1PkPj},i<jm[i, j]=\begin{cases} 0, & i=j\\ \min_{i\leqslant k<j}\{m[i, k]+m[k+1, j]+P_{i-1}P_kP_j\}, & i<j \end{cases}

该问题是满足最优子结构的。


利用数学归纳法,如果不做任何其他处理,可以证明上面动态规划方法的时间复杂度T(n)⩾2n−1T(n)\geqslant 2^{n-1},还是指数级别!

这是因为,每次遇到m[i,j]m[i, j],我们都当做一个子问题重新计算。实际上,可以发现我们进行了大量的重复计算,有很多子问题是重复的。所以,我们可以将子问题的解保存在一个表(备忘录)中,之后再遇到的时候,直接查表即可。

追踪解时,只需记录下每次最终选取的kk值即可。

记忆化之后,时间复杂度降低到O(n3)O(n^3)。

3.3、投资问题

有mm元钱,nn个投资项目,fi(x)f_i(x)是将xx元投入第ii个项目的收益。求出使得总收益最大的投资方案。


记Fk(x)F_k(x)是将xx元投入前kk个项目的最大收益。考虑将一部分前投入前k−1k-1个项目,剩下的前投给第kk个项目,则有递推式:

Fk(x)=max⁡0⩽y⩽x{Fk−1(x−y)+fk(y)}F_k(x)=\max_{0\leqslant y\leqslant x}\{F_{k-1}(x-y)+f_k(y)\}

时间复杂度为O(nm2)O(nm^2)。

3.4、一般背包问题

假设将nn种物品(每种物品有无数个)放入一个背包,第ii个物品的重量为wiw_i,价值为viv_i,背包的重量限制为bb。求出使得背包中物品的总价值最大的方案。


记只考虑前kk种物品,总重不超过yy时的最大价值为Fk(y)F_k(y)。每次考虑装至少装入一个物品kk还是不装,容易写出递推式:

Fk(y)=max⁡{Fk−1(y),Fk(y−wk)+vk}F_k(y)=\max\{F_{k-1}(y), F_{k}(y-w_k)+v_k\}

追踪解,只需同时记录一个ik(y)i_k(y),表示计算Fk(y)F_k(y)时,最终方案中装入物品的最大标号。具体来说,其更新公式如下:

ik(y)={ik−1(y),Fk−1(y)>Fk(y−wk)+vkk,Fk−1(y)⩽Fk(y−wk)+vki_k(y)=\begin{cases} i_{k-1}(y), & F_{k-1}(y)> F_k(y-w_k)+v_k\\ k, & F_{k-1}(y)\leqslant F_k(y-w_k)+v_k \end{cases}

Fn(b)F_n(b)即为最终的解。x=in(b)x=i_n(b)是最后一个装入的物品,in(b−wx)i_n(b-w_x)是倒数第二个装入的物品,以此类推。

时间复杂度为O(nb)O(nb)。

3.5、最长公共子序列问题

设两个序列X={x1,x2,…,xm}X=\{x_1, x_2, \dots, x_m\}和Z={z1,z2,…,zn}Z=\{z_1, z_2, \dots, z_n\}。如果存在j1<j2<⋯<jnj_1<j_2<\dots<j_n,使得zi=xji,∀i=1,2,…,nz_i=x_{j_i}, \forall i=1, 2, \dots, n,则称ZZ是XX的一个子序列。如果ZZ同时是XX和YY的子序列,就称它是XX和YY的公共子序列(LCS)。

给定两个序列XX和YY,求出它们的最长公共子序列。


蛮力算法。依次检查XX的所有子序列是否在YY中。子序列有2∣X∣2^{|X|}种,检查一个子序列是否存在另一个序列中,需要∣Y∣|Y|的时间。假定m=∣X∣⩽∣Y∣=nm=|X|\leqslant |Y|=n,则时间复杂度可写为O(n2m)O(n 2^m)


动态规划算法。记考虑XX的前ii个元素,以及YY的前jj个元素时的 LCS 长度为C[i,j]C[i, j]。可以想到,如果xi=yjx_i=y_j,则可以添加到目前 LCS 的末尾,使 LCS 长度增加。所以有递推公式如下:

C[i,j]={C[i−1,j−1]+1,xi=yjmax⁡{C[i−1,j],C[i,j−1]},xi≠xjC[i, j]=\begin{cases} C[i-1, j-1]+1, & x_i=y_j\\ \max\{C[i-1, j], C[i, j-1]\}, & x_i\neq x_j \end{cases}

追踪解,根据上面三种情况反向追踪即可。

总的时间复杂度为O(mn)O(mn)。

3.6、黑白图像存储问题

设图像像素序列为{p1,p2,…,pn}\{p_1, p_2, \dots, p_n\},其中每个像素点pip_i都是一个 0~255 的灰度值,需要 8 位来存储。

考虑对像素点序列进行分段S1,S2,…,SmS_1, S_2, \dots, S_m,段SiS_i有l[i]l[i]个像素(l[i]−1⩽255l[i]-1\leqslant 255),每个像素都用b[i]b[i]位来存储,则总的存储位数为(l[i]b[i]+8+3)⋅m(l[i]b[i]+8+3)\cdot m。其中bib_i满足:

b[i]=⌈log⁡2(max⁡pj∈Sipj+1)⌉b[i]=\left\lceil\log_2\left(\max_{p_j\in S_i}p_j+1\right)\right\rceil

求出最佳分段方案,使得总的存储位数最少。


动态规划算法。记S[i]S[i]是前ii个像素的采用最佳分段方案所需的存储位数,考虑最后一个段SmS_m如何划分,则有:

S[i]=min⁡1⩽j⩽min⁡{i,256}{S[i−j]+j⋅b[m]+11},bm=⌈log⁡2(max⁡pj∈Smpj+1)⌉S[i]=\min_{1\leqslant j\leqslant \min\{i, 256\}}\{S[i-j]+j\cdot b[m]+11\}, b_m=\left\lceil\log_2\left(\max_{p_j\in S_m}p_j+1\right)\right\rceil

追踪解只需要记录每次最终选择的jj值,这就是每段的长度。

总的时间复杂度为O(n)O(n)。

3.7、最大子串和问题

给定nn个数的序列(可能存在负数)A={a1,a2,…,an}A=\{a_1, a_2, \dots, a_n\},求出一个连续子串,使得子串的和最大。


蛮力算法。暴力枚举每一个子串,使用前缀和数组的情况下,时间复杂度为O(n2)O(n^2)。


分治算法。考虑前半段中的最大子串和以及后半段的最大子串和,以及跨越两个段的最大子串和。跨越两个段的情况,可以从中间往两边拓展,可在O(n)O(n)时间内求解。

所以时间复杂度T(n)=2T(n/2)+O(n)T(n)=2T(n/2)+O(n),由主定理,有T(n)=O(nlog⁡n)T(n)=O(n\log n)。


动态规划算法。记S[i]S[i]是以aia_i结尾的最大子串和,容易得到递推公式:

S[i]=max⁡{S[i−1]+ai,ai}S[i]=\max\{S[i-1]+a_i, a_i\}

时间复杂度为O(n)O(n)。

3.8、最优二叉搜索树问题

假设有形如下图的二叉搜索树:

BST

其中圆形节点表示是实际数据,方形节点是虚拟节点,是不在二叉搜索树中的数据最终落在的位置。每个节点都有一个概率,表示是搜索中最终节点的概率。

设实际数据集为{x1,x2,…,xn}\{x_1, x_2, \dots, x_n\},求出一种最优的二叉搜索树,使得搜索的期望比较次数最小。


动态规划算法。考虑某段数据{xi,xi+1,…,xj}\{x_i, x_{i+1}, \dots, x_j\},从中选择xkx_k作为根。记m[i,j]m[i, j]是这段数据的最优二叉搜索树的期望比较次数,w[i,j]w[i, j]为这段数据(包括实际数据以及虚拟节点)的概率和。提出根节点之后,到左右子树的比较次数都增加了 1。终点落在根节点上,需要一次比较。所以有:

m[i,j]=min⁡i⩽k⩽j{m[i,k−1]+1⋅w[i,k−1]+m[k+1,j]+1⋅w[k+1,j]+1⋅pk}m[i, j]=\min_{i\leqslant k\leqslant j}\{m[i, k-1]+1\cdot w[i, k-1] + m[k+1, j]+1\cdot w[k+1, j] + 1\cdot p_k\}

化简后,有:

m[i,j]=min⁡i⩽k⩽j{m[i,k−1]+m[k+1,j]+w[i,j]}m[i, j]=\min_{i\leqslant k\leqslant j}\{m[i, k-1]+m[k+1, j]+w[i, j]\}

总的时间复杂度为O(n3)O(n^3)。

第四章 贪心算法

1、基本概念

贪心法的基本思想是:在对问题求解时,总是做出在当前看来是最好的选择。也即,不从整体最优上加以考虑,只做局部最优解。显然,这样做并不一定能得到全局最优解。

贪心选择性质是指,一个问题的整体最优解可以通过一系列局部最优的选择得到。要想使用贪心算法得到最优解,必须证明问题具有贪心选择性质。

之前已经提到了最优子结构,一个问题拥有最优子结构是能够用动态规划算法以及贪心算法求解的关键特征。并不是所有具有最优子结构的问题都能够使用贪心算法求解,但是往往可以利用其来证明贪心选择性质。

2、数学归纳法

证明贪心选择性质时,常常用到数学归纳法。数学归纳法适合证明涉及自然数的命题P(n)P(n)。

2.1、第一数学归纳法

归纳基础:证明P(1)P(1)成立(或者P(0)P(0)成立)。
归纳步骤:假设对任意自然数kk,P(k)P(k)成立,证明P(k+1)P(k+1)成立。

2.2、第二数学归纳法

归纳基础:证明P(1),P(2),…,P(m)P(1), P(2), \dots, P(m)成立。
归纳步骤:假设对任意自然数kk,P(1),P(2),…,P(k)P(1), P(2), \dots, P(k)成立,证明P(k+1)P(k+1)成立。

3、实例

3.1、活动选择问题

设有nn个活动,活动ii的开始时间和结束时间分别为sis_i和fif_i。如果活动ii和活动jj满足si⩾fjs_i \geqslant f_j或者sj⩾fis_j \geqslant f_i,则称活动ii和活动jj是相容的,求出最大的两两相容的活动集合。


贪心策略 1:总是选择开始最早的活动。这样的话,如果某个活动持续非常久,会挤占掉其他活动,可能不如选稍晚一些开始,但很快结束的活动。

贪心策略 2:总是选持续时间最短的活动。这样的话,如果持续时间最短的活动开始的很晚,可能会错过很多活动。

贪心策略 3:总是选择结束时间最早的活动。先按照结束时间排序,然后扫描一遍,选出相容的活动。时间主要消耗在排序上,复杂度为O(nlog⁡n)O(n\log n)。


接下来证明贪心选择性质,下面的贪心策略代指贪心策略 3。假设活动集已按照结束时间升序排列,下面提及的序号都是排序后的序号。

命题:按照贪心策略已经选择了kk项活动i1=1,i2,…,iki_1=1, i_2, \dots, i_k,存在某个最优解AA包含活动i1=1,i2,…,iki_1=1, i_2, \dots, i_k。

归纳基础:证明k=1k=1时成立,也即证明最优解包含i1=1i_1=1。反证法,假设最优解中不包含i1=1i_1=1。任取一个最优解AA,将其中活动也按结束时间升序排列。由于活动 1 是结束时间最早的活动,一定有f1⩽fjf_1\leqslant f_j。用活动 1 替换掉活动jj,得到一个新的解A′A',对后续活动没有影响,也是一个最优解。这与假设矛盾,所以最优解中一定包含活动 1。

归纳步骤:假设对于kk,原命题成立,也即此时最优解AA包含i1=1,i2,…,iki_1=1, i_2, \dots, i_k。AA的剩余部分BB来自于集合S′={i∣i∈S,si⩾fik}S'=\{i| i\in S, s_i \geqslant f_{i_k}\}。由反证法易得,BB一定是S′S'的最优解。考虑活动集为S′S'的新问题,由归纳假设,S′S'的B′B'一定包含S′S'中结束时间最早的活动。由{i1,i2,…,ik}∪B′={i1,i2,…,ik,ik+1}∪(B′−{ik+1)})\{i_1, i_2, \dots, i_k\}\cup B'=\{i_1, i_2, \dots, i_k, i_{k+1}\} \cup (B'-\{i_{k+1})\})是最优解,可得最优解包含i1,i2,…,ik,ik+1i_1, i_2, \dots, i_k, i_{k+1}。所以原命题对于k+1k+1也成立。

第六章 线性规划

1、基本概念

1.1、一般形式

线性规划问题的一般形式如下:

min⁡(max⁡)z=∑j=1ncjxjs.t.∑j=1naijxj⩽(=,⩾)bi,i=1,2,…,mxj⩾0,j∈J⊆{1,2,…,n}xj,j∈{1,2,…,n}−J\begin{align*} \min(\max) z=\sum_{j=1}^n c_jx_j\\ \text{s.t.} \sum_{j=1}^n a_{ij}x_j \leqslant(=, \geqslant)b_i, & i=1, 2, \dots, m\\ x_j\geqslant 0,& j\in J\subseteq \{1, 2, \dots, n\}\\ x_j, &j\in \{1, 2, \dots, n\}-J \end{align*}

这四行从上到下依次代表:目标函数,约束条件,非负约束条件,自由变量。


可行解:满足约束条件和非负条件的变量
可行域:所有可行解的集合
最优解:在可行域中,目标函数取得最小(或最大)值的解
最优值:最优解对应的目标函数值

在高中已经学习过了二维线性规划的图解法,可行域由多条代表约束条件的直线围城,是一个凸多边形(可能无界,也可能是空集)。如果有最优解,一定在凸多边形的顶点取到。解的情况有四种:

  • 有唯一最优解
  • 有无穷多个最优解
  • 有可行解,但无最优解
  • 无可行解,也无最优解

推广到一般的nn维线性规划也是如此。

1.2、标准形

一般线性规划问题,总是可以写成标准形:

min⁡z=∑j=1ncjxjs.t.∑j=1naijxj=bi,i=1,2,…,mxj⩾0,j=1,2,…,n\begin{align*} \min z=\sum_{j=1}^n c_jx_j\\ \text{s.t.} \sum_{j=1}^n a_{ij}x_j=b_i, & i=1, 2, \dots, m\\ x_j\geqslant 0,& j=1, 2, \dots, n \end{align*}

将max⁡\max改min⁡\min,或者不等号变号比较简单,不再赘述。对于不等式约束条件,可以通过引入松弛变量/剩余变量,将其转化为等式约束条件:

∑j=1naijxj⩽bi⇒∑j=1naijxj+yi=bi,yi⩾0∑j=1naijxj⩾bi⇒∑j=1naijxj−yi=bi,yi⩾0\begin{align*} \sum_{j=1}^n a_{ij}x_j\leqslant b_i \Rightarrow \sum_{j=1}^n a_{ij}x_j+y_i=b_i, y_i\geqslant 0 \\ \sum_{j=1}^n a_{ij}x_j\geqslant b_i \Rightarrow \sum_{j=1}^n a_{ij}x_j-y_i=b_i, y_i\geqslant 0 \end{align*}

其中第一行引入的变量称为松弛变量,第二行引入的变量称为剩余变量。

对于自由变量:

xj∈R⇒xj=xj′−xj′′,xj′⩾0,xj′′⩾0x_j\in R\Rightarrow x_j=x_j'-x_j'', x_j'\geqslant 0, x_j''\geqslant 0

1.3、矩阵形式

标准形可以写为矩阵形式。其中目标函数可写为:

min⁡z=cTx=[c1c2⋮cn]T[x1x2⋮xn]\min z=\boldsymbol{c}^T\boldsymbol{x}=\begin{bmatrix} c_1 \\ c_2 \\ \vdots \\ c_n \end{bmatrix}^T \begin{bmatrix} x_1\\ x_2\\ \vdots\\ x_n \end{bmatrix}

约束条件可写为:

[a11a12…a1na21a22…a2n⋮⋮⋱⋮am1am2…amn][x1x2⋮xn]=[P1P2…Pn][x1x2⋮xn]=Ax=b\begin{bmatrix} a_{11} & a_{12} & \dots & a_{1n} \\ a_{21} & a_{22} & \dots & a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{m1} & a_{m2} & \dots & a_{mn} \end{bmatrix} \begin{bmatrix} x_1\\ x_2\\ \vdots\\ x_n \end{bmatrix}=\begin{bmatrix} \boldsymbol{P}_1 & \boldsymbol{P}_2 & \dots & \boldsymbol{P}_n \end{bmatrix} \begin{bmatrix} x_1\\ x_2\\ \vdots\\ x_n \end{bmatrix}=\boldsymbol{A}\boldsymbol{x}=\boldsymbol{b}

2、标准形的解

2.1、一些定义

设A\boldsymbol{A}的秩为rr,A\boldsymbol{A}的mm个线性无关的列向量称为标准型的基。

给定标准形的一组基B={Pi1,Pi2,…,Pim}\boldsymbol{B}=\{\boldsymbol{P}_{i_1}, \boldsymbol{P}_{i_2}, \dots, \boldsymbol{P}_{i_m}\},对应变量xi1,xi2,…,ximx_{i_1}, x_{i_2}, \dots, x_{i_m}称为基变量,其余变量称为非基变量。

基变量构成的向量记为xB\boldsymbol{x}_B,非基变量构成的向量记为xN\boldsymbol{x}_N。令xN=0\boldsymbol{x}_N=\boldsymbol{0},则等式约束变为BxB=b\boldsymbol{B}\boldsymbol{x}_B=\boldsymbol{b},解得xB=B−1b\boldsymbol{x}_B=\boldsymbol{B}^{-1}\boldsymbol{b}。将xB\boldsymbol{x}_B和xN\boldsymbol{x}_N重新组装为x\boldsymbol{x},这个x\boldsymbol{x}显然满足等式约束,且非基变量全为 0,称其是关于基B\boldsymbol{B}的基本解(系数矩阵满秩,所以基本解是唯一的)。如果x\boldsymbol{x}是基本解,且同时还满足非负约束xi⩾0,∀ix_i\geqslant 0, \forall i,则称其为基本可行解,对应的基称为一个可行基。

2.2、基本可行解的性质

引理:Ax=b\boldsymbol{Ax}=\boldsymbol{b}的解α\boldsymbol{\alpha}是基本解当且仅当α\boldsymbol{\alpha}的非零分量对应的列向量线性无关。

必要性,由基本解的定义立得。

充分性,设α\boldsymbol{\alpha}非零分量对应的列向量为Pj1,Pj2,…,Pjr\boldsymbol{P}_{j_1}, \boldsymbol{P}_{j_2}, \dots, \boldsymbol{P}_{j_r},它们是线性无关的。由于A\boldsymbol{A}的秩为mm,必然存在其他m−rm-r个列向量Pjr+1,Pjr+2,…,Pjm\boldsymbol{P}_{j_{r+1}}, \boldsymbol{P}_{j_{r+2}}, \dots, \boldsymbol{P}_{j_m},这一共mm个列向量线性无关,是A\boldsymbol{A}的基B\boldsymbol{B}。则α\boldsymbol{\alpha}也是方程BxB=b\boldsymbol{B}\boldsymbol{x}_B=\boldsymbol{b}的解,由解的唯一性,α\boldsymbol{\alpha}就是基本解。


定理 1:若标准形有可行解,则必有基本可行解。

证明:设α\boldsymbol{\alpha}是一个可行解。设其非零分量为α1,α2,…,αr\alpha_1, \alpha_2, \dots, \alpha_r,对应的列向量为Pj1,Pj2,…,Pjr\boldsymbol{P}_{j_1}, \boldsymbol{P}_{j_2}, \dots, \boldsymbol{P}_{j_r}。由引理,若这rr个列向量线性无关,所以α\boldsymbol{\alpha}是基本解。

若不然,由线性无关,存在不全为 0 的λ1,λ2,…,λr\lambda_1, \lambda_2, \dots, \lambda_r,使得∑i=1rλiPji=0\sum_{i=1}^r \lambda_i\boldsymbol{P}_{j_i}=\boldsymbol{0}。取λr+1=λr+2=⋯=λn=0\lambda_{r+1}=\lambda_{r+2}=\dots=\lambda_n=0,则有∑i=1nλiPji=0\sum_{i=1}^n \lambda_i\boldsymbol{P}_{j_i}=\boldsymbol{0}。于是,对任意的δ\delta,有:

∑i=1n(αi+δλi)Pji=∑i=1nαiPji+δ∑i=1nλiPji=b\sum_{i=1}^n (\alpha_i+\delta \lambda_i)\boldsymbol{P}_{j_i}=\sum_{i=1}^n \alpha_i\boldsymbol{P}_{j_i}+\delta\sum_{i=1}^n \lambda_i\boldsymbol{P}_{j_i}=\boldsymbol{b}

如果想让α+δλ\boldsymbol{\alpha}+\delta \boldsymbol{\lambda}成为一个可行解,则αi+δλi⩾0,∀i\alpha_i+\delta \lambda_i\geqslant 0, \forall i

  • λi=0\lambda_i=0时,恒成立
  • λi>0\lambda_i>0时,则要有δ⩾−αiλi\delta \geqslant -\frac{\alpha_i}{\lambda_i};当λj<0\lambda_j<0时,要有δ⩽−αiλi\delta \leqslant -\frac{\alpha_i}{\lambda_i}。设k=arg⁡min⁡i,λi≠0{∣αiλi∣}k=\arg\min_{i, \lambda_i\neq 0}\{|\frac{\alpha_i}{\lambda_i}|\},则取δ∗=−αkλk\delta^*= -\frac{\alpha_k}{\lambda_k}。显然β=α+δ∗λ\boldsymbol{\beta}=\boldsymbol{\alpha}+\delta^*\boldsymbol{\lambda}也是一个可行解,而且比α\boldsymbol{\alpha}少一个非零分量(αk+δ∗λk=0\alpha_k+\delta^*\lambda_k=0)。重复上述过程至多rr次,就得到了一个基本可行解。

综上,证毕。


定理 2:若标准形有最优解,则必定存在一个基本可行解是最优解。

证明:只需在定理 1 的基础上,证明α\boldsymbol{\alpha}是最优解时,β\boldsymbol{\beta}也是最优解。

设α\boldsymbol{\alpha},显然其也是可行解。对于其任意零分量αi\alpha_i,一定有λi=0\lambda_i=0,所以对于任意δ\delta,有αi±δλi⩾0\alpha_i\pm \delta \lambda_i\geqslant 0;对于任意非零分量αi\alpha_i,有λi≠0\lambda_i\neq 0,取一个足够小的δ>0\delta>0,使仍然有αi±δλi⩾0\alpha_i\pm \delta \lambda_i\geqslant 0,且等式约束也满足。所以α±δλ\boldsymbol{\alpha}\pm\delta \boldsymbol{\lambda}也是一个可行解。由α\boldsymbol{\alpha}是最优解,于是有:

∑i=1nciαi⩽∑i=1nci(αi±δλi)=∑i=1nciαi±δ∑i=1nciλi\sum_{i=1}^n c_i\alpha_i\leqslant \sum_{i=1}^n c_i(\alpha_i\pm \delta \lambda_i)=\sum_{i=1}^n c_i\alpha_i\pm \delta \sum_{i=1}^n c_i\lambda_i

可得∑i=1nciλi=0\sum_{i=1}^n c_i\lambda_i=0,所以:

∑i=1nciβi=∑i=1nci(αi+δ∗λi)=∑i=1nciαi\sum_{i=1}^n c_i\beta_i=\sum_{i=1}^n c_i(\alpha_i+ \delta^* \lambda_i)=\sum_{i=1}^n c_i\alpha_i

所以β\boldsymbol{\beta}也是最优解。所以按照定理 1 中的流程,最终推得的基本可行解也是最优解,定理 2 得证。


综上,要找到原问题的一个最优解,在标准形中的基本可行解中寻找即可。A\boldsymbol{A}至多有CnmC_{n}^{m}个基,故至多有CnmC_{n}^{m}个基本可行解,这就是我们的搜索空间。

3、单纯形法

基本步骤如下:

  1. 选取一个初始可行基,确定初始基本可行解
  2. 检查当前的基本可行解。若是最优解或发现无最优解,则结束;否则作基变换,用一个非基变量替换一个基变量,得到新的基和对应的基本可行解,且使目标函数值至少不增。

3.1、确定初始基本可行解

考虑最简单的情况,设约束条件为:

∑j=1naijxj⩽0,i=1,2,…,m\sum_{j=1}^n a_{ij}x_j \leqslant 0, i=1, 2, \dots, m

其中bi⩾0b_i\geqslant 0。考虑引入mm个松弛变量xn+i⩾0x_{n+i}\geqslant 0,约束变为:

∑j=1naijxj+xn+i=0,i=1,2,…,m\sum_{j=1}^n a_{ij}x_j +x_{n+i}=0, i=1, 2, \dots, m

选取{xn+i}\{x_{n+i}\}作为基向量,易得其基本可行解为(0,0,…,0,b1,b2,…,bm)T(0, 0, \dots, 0, b_1, b_2, \dots, b_m)^T

3.2、最优性检验

考虑某个可行基B=(Pπ(1),Pπ(2),⋯ ,Pπ(m))\boldsymbol{B}=(\boldsymbol{P}_{\pi (1)}, \boldsymbol{P}_{\pi (2)}, \cdots, \boldsymbol{P}_{\pi (m)})(这里的π(⋅)\pi(\cdot)是一种映射,表示B\boldsymbol{B}的第ii个列向量对应着A\boldsymbol{A}的第π(i)\pi(i)个列向量)。记A\boldsymbol{A}中非基变量的列构成的矩阵为N\boldsymbol{N},有:

B−1Ax=B−1(BxB+NxN)=xB+B−1NxN=B−1b\boldsymbol{B}^{-1}\boldsymbol{A}\boldsymbol{x}=\boldsymbol{B}^{-1}(\boldsymbol{B}\boldsymbol{x}_B+\boldsymbol{N}\boldsymbol{x}_N)=\boldsymbol{x}_B+\boldsymbol{B}^{-1}\boldsymbol{N}\boldsymbol{x}_N=\boldsymbol{B}^{-1}\boldsymbol{b}

可解得xB=B−1b−B−1NxN\boldsymbol{x}_B=\boldsymbol{B}^{-1}\boldsymbol{b}-\boldsymbol{B}^{-1}\boldsymbol{N}\boldsymbol{x}_N

代入目标函数:

z=cTx=cBTxB+cNTxN=cBTB−1b+(cNT−cBTB−1N)xNz=\boldsymbol{c}^T\boldsymbol{x}=\boldsymbol{c}_B^T\boldsymbol{x}_B+\boldsymbol{c}_N^T\boldsymbol{x}_N=\boldsymbol{c}_B^T\boldsymbol{B}^{-1}\boldsymbol{b}+(\boldsymbol{c}_N^T-\boldsymbol{c}_B^T\boldsymbol{B}^{-1}\boldsymbol{N})\boldsymbol{x}_N

所以,对于基B\boldsymbol{B},其基本可行解x\boldsymbol{x}由xB=B−1b\boldsymbol{x}_B=\boldsymbol{B}^{-1}\boldsymbol{b}以及xN=0\boldsymbol{x}_N=\boldsymbol{0}构成,对应的目标函数值为z0=cBTB−1bz_0=\boldsymbol{c}_B^T\boldsymbol{B}^{-1}\boldsymbol{b}


将z0z_0代入,继续运算,得到简化的目标函数:

z=z0+(cNT−cBTB−1N)xN=z0+(cBT−cBTB−1B)xB+(cNT−cBTB−1N)xN=z0+(cT−cBTB−1A)x=z0+λTx\begin{align*} z&=z_0+(\boldsymbol{c}_N^T-\boldsymbol{c}_B^T\boldsymbol{B}^{-1}\boldsymbol{N})\boldsymbol{x}_N=z_0+(\boldsymbol{c}_B^T-\boldsymbol{c}_B^T\boldsymbol{B}^{-1}\boldsymbol{B})\boldsymbol{x}_B+(\boldsymbol{c}_N^T-\boldsymbol{c}_B^T\boldsymbol{B}^{-1}\boldsymbol{N})\boldsymbol{x}_N \\ &=z_0+(\boldsymbol{c}^T-\boldsymbol{c}_B^T\boldsymbol{B}^{-1}\boldsymbol{A})\boldsymbol{x}=z_0+\boldsymbol{\lambda}^T\boldsymbol{x} \end{align*}

其中λT=cT−cBTB−1A\boldsymbol{\lambda}^T=\boldsymbol{c}^T-\boldsymbol{c}_B^T\boldsymbol{B}^{-1}\boldsymbol{A}。

称λ\boldsymbol{\lambda}的分量λ1,λ2,…,λn\lambda_1, \lambda_2, \dots, \lambda_n为检验数,对应基变量的检验数必为 0。


记B−1A=α\boldsymbol{B}^{-1}\boldsymbol{A}=\boldsymbol{\alpha},β=B−1b\boldsymbol{\beta}= \boldsymbol{B}^{-1}\boldsymbol{b},Pj′=B−1Pj=(α1j,α2j,…,αmj)T\boldsymbol{P}_j'=\boldsymbol{B}^{-1}\boldsymbol{P}_j=(\alpha_{1j}, \alpha_{2j}, \dots, \alpha_{mj})^T

定理 3:对于可行基B\boldsymbol{B},给定其基本可行解x(0)\boldsymbol{x}^{(0)},若λi⩾0,∀i\lambda_i\geqslant 0, \forall i,则x(0)\boldsymbol{x}^{(0)}是最优解;若存在λk<0\lambda_k<0且所有αik⩽0\alpha_{ik}\leqslant 0,则原问题无最优解。

证明:如果λi⩾0,∀i\lambda_i \geqslant 0, \forall i,则对任意可行解x\boldsymbol{x},必有λTx⩾0\boldsymbol{\lambda}^T\boldsymbol{x}\geqslant 0,所以z⩾z0z\geqslant z_0,x(0)\boldsymbol{x}^{(0)}是最优解。

如果存在λk<0\lambda_k<0,由λ\boldsymbol{\lambda}的定义,λk\lambda_k一定对应某个非基变量。取xk=M>0x_k=M>0,其他非基变量取 0,可求得xBi=βi−αikM⩾0x_{Bi}=\beta_i-\alpha_{ik}M\geqslant 0,所以这也是一个可行解。其对应的目标函数值为z=z0+λkM+Cz=z_0+\lambda_k M+C,其中CC是基变量对应的分量,有C⩾0C\geqslant 0。当M→+∞M\to +\infty时,z→−∞z\to -\infty,所以原问题无最优解。


还有一种可能的情况是,存在λk<0\lambda_k<0,但是αik\alpha_{ik}不全为非正数。这种情况,就要用到下面介绍的基变换,来进一步进行下去。

3.3、基变换

设λk<0\lambda_k<0且存在αlk>0\alpha_{lk}>0,其对应的xkx_k一定是非基变量。进行基变换:用非基变量xkx_k替换基变量xπ(l)x_{\pi (l)},对应得到新的基B′={Pπ(1),…,Pπ(l−1),Pk,Pπ(l+1),…,Pπ(m)}\boldsymbol{B}'=\{\boldsymbol{P}_{\pi (1)}, \dots, \boldsymbol{P}_{\pi (l-1)}, \boldsymbol{P}_k, \boldsymbol{P}_{\pi (l+1)}, \dots, \boldsymbol{P}_{\pi (m)}\}。称xk\boldsymbol{x}_k为换入变量,xπ(l)\boldsymbol{x}_{\pi (l)}为换出变量。

首先要证明B′\boldsymbol{B}'确实是一个基。只需证明Pπ(l)\boldsymbol{P}_{\pi (l)}可以被B′\boldsymbol{B}'表示即可(因为B\boldsymbol{B}是一组基,将这组基下的线性表示中的Pπ(l)\boldsymbol{P}_{\pi (l)}换为B′\boldsymbol{B}'下的表示,就得到了纯B′\boldsymbol{B}'下的线性表示)。

由于(Pπ(1)′,Pπ(2)′,…,Pπ(m)′)=B−1B=I(\boldsymbol{P}_{\pi (1)}', \boldsymbol{P}_{\pi (2)}', \dots, \boldsymbol{P}_{\pi (m)}')=\boldsymbol{B}^{-1}\boldsymbol{B}=\boldsymbol{I},所以有Pk′=∑i=1mαikPπ(i)′\boldsymbol{P}_k'=\sum_{i=1}^m \alpha_{ik}\boldsymbol{P}_{\pi (i)}'。在两边同左乘B\boldsymbol{B},移项化简,最终可得:

Pπ(l)=1αlkPk−∑i≠lαikαlkPπ(i)\boldsymbol{P}_{\pi (l)}=\frac{1}{\alpha_{lk}}\boldsymbol{P}_k-\sum_{i\neq l}\frac{\alpha_{ik}}{\alpha_{lk}}\boldsymbol{P}_{\pi (i)}

因此,B′\boldsymbol{B}'确实是一个基。


由(Pπ(1)′,Pπ(2)′,…,Pπ(m)′)=B−1B=I(\boldsymbol{P}_{\pi (1)}', \boldsymbol{P}_{\pi (2)}', \dots, \boldsymbol{P}_{\pi (m)}')=\boldsymbol{B}^{-1}\boldsymbol{B}=\boldsymbol{I},将单位阵的第ii列换成Pk′\boldsymbol{P}_k'得到H\boldsymbol{H},再左乘B\boldsymbol{B},实际上就得到了B′\boldsymbol{B}'对应的矩阵。也即有B′=BH\boldsymbol{B}'=\boldsymbol{B}\boldsymbol{H}。于是有:

B′−1Ax=B′−1b⇔H−1B−1Ax=H−1B−1b=H−1β\boldsymbol{B}'^{-1}\boldsymbol{A}\boldsymbol{x}=\boldsymbol{B}'^{-1}b\Leftrightarrow \boldsymbol{H}^{-1}\boldsymbol{B}^{-1}\boldsymbol{A}\boldsymbol{x}=\boldsymbol{H}^{-1}\boldsymbol{B}^{-1}b=\boldsymbol{H}^{-1}\boldsymbol{\beta}

所以用B′\boldsymbol{B}'代替B\boldsymbol{B},等价于在原来的基础上左乘H−1\boldsymbol{H}^{-1}。H−1\boldsymbol{H}^{-1}具体如下:

H−1=[1  −α1k/αlk    ⋱ ⋮    1−αl−1,k/αlk     1/αlk     −αl+1,k/αlk1    ⋮ ⋱   −αmk/αlk  1]\boldsymbol{H}^{-1}=\begin{bmatrix} 1 & \ & \ & -\alpha_{1k}/\alpha_{lk} & \ & \ & \ \\ \ & \ddots & \ &\vdots & \ & \ \\ \ & \ & 1 & -\alpha_{l-1,k}/\alpha_{lk} & \ & \ \\ \ & \ & \ & 1/\alpha_{lk} & \ & \ \\ \ & \ & \ & -\alpha_{l+1,k}/\alpha_{lk} & 1 & \ \\ \ & \ & \ & \vdots & \ & \ddots \\ \ & \ & \ & -\alpha_{mk}/\alpha_{lk} & \ & \ & 1 \end{bmatrix}


记H−1B−1A=[αij′]m×n\boldsymbol{H}^{-1}\boldsymbol{B}^{-1}\boldsymbol{A}=\begin{bmatrix} \alpha'_{ij} \end{bmatrix}_{m\times n},β′=H−1B−1b\boldsymbol{\beta}'= \boldsymbol{H}^{-1}\boldsymbol{B}^{-1}\boldsymbol{b}

于是可以写出αij′\alpha'_{ij}和βi′\beta_i'的表达式:

αij′={αljαlk,i=lαij−αikαljαlk,i≠l\alpha'_{ij}=\begin{cases} \frac{\alpha_{lj}}{\alpha_{lk}} , & i=l \\ \alpha_{ij}-\frac{\alpha_{ik}\alpha_{lj}}{\alpha_{lk}} , & i\neq l \end{cases}

βi′={βlαlk,i=lβi−αikβlαlk,i≠l\beta_i'=\begin{cases} \frac{\beta_l}{\alpha_{lk}} , & i=l \\ \beta_i-\frac{\alpha_{ik}\beta_l}{\alpha_{lk}} , & i\neq l \end{cases}

直观理解,就是把变形后的约束方程xB+B−1NxN=B−1b\boldsymbol{x}_B+\boldsymbol{B}^{-1}\boldsymbol{N}\boldsymbol{x}_N=\boldsymbol{B}^{-1}\boldsymbol{b}中第ll个方程中xlx_l的系数变为 1,然后用这个方程消去其他方程中的xlx_l。

要保证B′\boldsymbol{B}'是可行的,需要保证βi′⩾0,∀i\beta_i'\geqslant 0, \forall i。由βi⩾0,∀i\beta_i\geqslant 0, \forall i以及αlk>0\alpha_{lk}>0,所以αik⩽0\alpha_{ik}\leqslant 0时,有βi′⩾0\beta_i'\geqslant 0成立。而αlk>0\alpha_{lk}>0,要有:

βlαlk⩽βiαik\frac{\beta_l}{\alpha_{lk}}\leqslant \frac{\beta_i}{\alpha_{ik}}

因此,只需取l=arg⁡min⁡i,αik>0{βi/αik}l=\arg\min_{i, \alpha_{ik}>0}\{ \beta_i/\alpha_{ik}\}


相应的,对化简的目标函数也做类似的变换,得到基于B′\boldsymbol{B}'的简化目标函数,其中:

λj′=λj−λkαljαlkz0′=z0+λkβlαlk\begin{align*} \lambda_j'&=\lambda_j-\frac{\lambda_k\alpha_{lj}}{\alpha_{lk}}\\ z_0'&=z_0+\frac{\lambda_k\beta_l}{\alpha_{lk}} \end{align*}

直观上理解,用之前提到的化简后的第ll个方程,消去z=z0+λTxz=z_0+\boldsymbol{\lambda}^T\boldsymbol{x}中的xkx_k。

3.4、单纯形法的完整叙述

  1. 设初始可行基B=(Pπ(1),Pπ(2),⋯ ,Pπ(m))\boldsymbol{B}=(\boldsymbol{P}_{\pi (1)}, \boldsymbol{P}_{\pi (2)}, \cdots, \boldsymbol{P}_{\pi (m)}),α=B−1A\boldsymbol{\alpha}= \boldsymbol{B}^{-1}\boldsymbol{A},β=B−1b\boldsymbol{\beta}= \boldsymbol{B}^{-1}\boldsymbol{b},λT=cT−cBTα\boldsymbol{\lambda}^T=\boldsymbol{c}^T-\boldsymbol{c}_B^T\boldsymbol{\alpha},z0=cBTβz_0=\boldsymbol{c}_B^T\boldsymbol{\beta}
  2. 若λj⩾0,1⩽j⩽n\lambda_j\geqslant 0, 1\leqslant j\leqslant n,则xB=β,xN=0\boldsymbol{x}_B=\boldsymbol{\beta}, \boldsymbol{x}_N=\boldsymbol{0}组合出的x\boldsymbol{x}是最优解,算法结束
  3. 否则,任取一个λk<0\lambda_k<0。若所有αik⩽0\alpha_{ik}\leqslant 0,则原问题无最优解,算法结束
  4. 否则存在ll使得αlk>0\alpha_{lk}>0,其中l=arg⁡min⁡i,αik>0{βi/αik}l=\arg\min_{i, \alpha_{ik}>0}\{ \beta_i/\alpha_{ik}\}
  5. 以xkx_k为换入变量,xπ(l)x_{\pi (l)}为换出变量,进行基变换。回到步骤 2

3.5、单纯形表

单纯形表如下图所示:

table

首先,来解释一下单纯形表的结构:

  • 最中心的大矩形存放矩阵α\boldsymbol{\alpha}
  • xB\boldsymbol{x}_B列存放基变量,cB\boldsymbol{c}_B存放的是基变量对应的目标函数系数,b\boldsymbol{b}列存放基变量的值
  • 最后一行存放检验数λ\boldsymbol{\lambda}
  • 将简化的目标函数改写为−z+∑j=1nλjxj=−z0-z+\sum_{j=1}^n \lambda_j \boldsymbol{x}_j=-z_0,与∑j=1nαijxj=βi\sum_{j=1}^n \alpha_{ij}x_j=\beta_i的形式统一起来。每一行就代表了其中一个等式。

由于选择的初始基B\boldsymbol{B}是单位阵,所以一开始α=A\boldsymbol{\alpha}=\boldsymbol{A},cB=0\boldsymbol{c}_B=\boldsymbol{0},λ=c\boldsymbol{\lambda}=\boldsymbol{c}


接下来讲述怎么利用单纯形表进行计算。

  1. 初始化。按照上面的叙述,初始化单纯形表

  2. 检查λj\lambda_j,如果全是非负,则已经得到最优解(由当前cB\boldsymbol{c}_B列和b\boldsymbol{b}列每行逐对相乘再求和得到,也就是z0z_0值)。如果存在λk<0\lambda_k<0(有多个时,选∣λk∣|\lambda_k|最大的那个),检查α\boldsymbol{\alpha}的第kk列,如果该列全为非正数,则无最优解。

  3. 否则,对每一个αik>0\alpha_{ik}>0,计算θi=βi/αik\theta_i=\beta_i/\alpha_{ik},填到θ\boldsymbol{\theta}列中。假设θl\theta_l最小,将αlk\alpha_{lk}圈起来,取xkx_k为换入变量,xπ(l)x_{\pi (l)}为换出变量,进行基变换。

  4. 在单纯形表上做基变换的操作:

    • 更新α,β,λ\boldsymbol{\alpha}, \boldsymbol{\beta}, \boldsymbol{\lambda}。做一些行变换,在 3.3 基变换章节已经叙述过了。
    • 将xB\boldsymbol{x}_B列的xπ(l)\boldsymbol{x}_{\pi (l)}换成xk\boldsymbol{x}_k。将cB\boldsymbol{c}_B的cπ(l)\boldsymbol{c}_{\pi (l)}换成ck\boldsymbol{c}_k。
  5. 这样,就得到了一张新的单纯形表,回到步骤 2

4、人工变量和两阶段法

4.1、人工变量

在 3.1 中,只讨论了对最简单的情况如何确定初始基本可行解。约束条件还有其他两种情况:

  • ∑j=1naijxj⩾bi\sum_{j=1}^n a_{ij}x_j \geqslant b_i
  • ∑j=1naijxj=bi\sum_{j=1}^n a_{ij}x_j = b_i

对于情况 1,引入剩余变量可以转化为情况 2。对于情况 2,引入人工变量yi⩾0y_i\geqslant 0,将约束变为:

∑j=1naijxj+yi=bi\sum_{j=1}^n a_{ij}x_j +y_i = b_i

一开始取所有的松弛变量和人工变量得到初始可行基,以及相应初始基本可行解。

4.2、两阶段法

引入人工变量后求解,最终得到的解中,只有所有的人工变量值均为 0 时,才能直接舍弃掉人工变量,得到原问题的解。否则,是不满足约束条件的。

设原问题为:

min⁡z=∑j=1ncjxjs.t.∑j=1naijxj=bi,1⩽i⩽mxj⩾0,1⩽j⩽n\begin{align*} \min z=\sum_{j=1}^n c_jx_j\\ \text{s.t.} \sum_{j=1}^n a_{ij}x_j = b_i, & 1\leqslant i\leqslant m\\ x_j\geqslant 0,& 1\leqslant j\leqslant n \end{align*}

引入人工变量xn+i,1⩽i⩽mx_{n+i}, 1\leqslant i\leqslant m,引入一个辅助问题:

min⁡w=∑i=1mxn+is.t.∑j=1naijxj+xn+i=bi,1⩽i⩽mxj⩾0,1⩽j⩽n+m\begin{align*} \min w=\sum_{i=1}^m x_{n+i}\\ \text{s.t.} \sum_{j=1}^n a_{ij}x_j +x_{n+i} = b_i, & 1\leqslant i\leqslant m\\ x_j\geqslant 0,& 1\leqslant j\leqslant n+m \end{align*}

由于w⩾0w\geqslant 0,所以辅助问题必有最优解。设其最优解为(x1∗,x2∗,…,xn+m∗)(x_1^*, x_2^*, \dots, x_{n+m}^*),最优值为w∗w^*,有以下三种情况:

  • w∗>0w^*>0,则原问题无可行解
  • w∗=0w^*=0且最优解中所有人工变量都是非基变量,也即xn+i∗=0,1⩽i⩽mx_{n+i}^*=0, 1\leqslant i\leqslant m,则(x1∗,x2∗,…,xn∗)(x_1^*, x_2^*, \dots, x_n^*)是原问题的基本可行解
  • w∗=0w^*=0但最优解中有人工变量是基变量。仍有xn+i∗=0,1⩽i⩽mx_{n+i}^*=0, 1\leqslant i\leqslant m,不妨先假设只有xn+kx_{n+k}是基变量。考虑解该辅助问题时最终的单纯形表中,xn+kx_{n+k}对应的行,假设在第rr行。由αx=β\boldsymbol{\alpha x}=\boldsymbol{\beta}有:

    ∑j=1m+nαrjxj=βr⇔∑j=1nαrjxj+1⋅xn+k+∑1⩽j⩽m,j≠kαr,n+jxn+j=0\sum_{j=1}^{m+n} \alpha_{rj}x_j=\beta_r \Leftrightarrow \sum_{j=1}^{n} \alpha_{rj}x_j+1\cdot x_{n+k}+ \sum_{1\leqslant j\leqslant m, j\neq k} \alpha_{r, {n+j}}x_{n+j}=0

    说明原问题中mm个约束等式并不是线性无关的(用单纯形表进行计算,本质上是不停的对[αβ]\begin{bmatrix} \boldsymbol{\alpha} & \boldsymbol{\beta} \end{bmatrix}进行行变换,这里组合出了非平凡的零解)。
    • 如果αrj,1⩽j⩽n\alpha_{rj}, 1\leqslant j\leqslant n全为 0,说明引入xn+kx_{n+k}的那个约束是冗余的,可以直接删掉
    • 如果存在αrl≠0\alpha_{rl} \neq 0,以xlx_{l}为换入变量,x∗n+kx*{n+k}为换出变量,进行基变换。由于β∗r=0\beta*{r}=0,所以这个基变换不会改变β\boldsymbol{\beta}列,也不会改变ww的值,但是使基变量中的人工变量变少了。

所以如果w∗=0w^*=0,最终总能变成情况 2,保证基变量中没有人工变量。

5、单纯形法的有限终止

如果基本可行解中基变量的值都大于 0, 则称这个基本可行解是非退化的, 否则称作退化的。

如果线性规划的所有基本可行解都是非退化的, 则称这个线性规划是非退化的。如果线性规划有可行解并且是非退化的, 则每一次基变换都会使目标函数值严格下降,从而在计算过程中,可行基不会从夫出现,因此单纯形法一定会在有限步内终止。

Bland 规则指出:

  • 当有多个λj<0\lambda_j<0时,去对应的非基变量中下标最小的那个作为换入变量
  • 当有多个θi=βi/αik\theta_i=\beta_i/\alpha_{ik}同时取到最小值时,取对应的基变量中下标最小的那个作为换出变量

按照 Bland 规则选取换入换出变量,可以保证单纯形法的有限终止。

6、对偶性

设线性规划:

max⁡cTxs.t.Ax⩽bx⩾0\begin{align*} \max \boldsymbol{c}^T\boldsymbol{x}\\ \text{s.t.} \boldsymbol{A}\boldsymbol{x}\leqslant \boldsymbol{b}\\ \boldsymbol{x}\geqslant \boldsymbol{0} \end{align*}

其对偶线性规划(也简称对偶/对偶规划)为:

min⁡bTys.t.ATy⩾cy⩾0\begin{align*} \min \boldsymbol{b}^T\boldsymbol{y}\\ \text{s.t.} \boldsymbol{A}^T\boldsymbol{y}\geqslant \boldsymbol{c}\\ \boldsymbol{y}\geqslant \boldsymbol{0} \end{align*}

第七章 网络流算法

1、有向图

图、节点、边的概念不在此赘述。

有向图可记作G=(V,E)G=(V, E),其中VV是节点集合,EE是边集合。每条边e∈Ee\in E可以表示为一个有序对(i,j)(i, j),其中i,j∈Vi, j\in V,表示从节点ii到节点jj的一条边。

从一个节点出发,到达另一个节点,所经过的边的序列称为一条路径。路径上边的个数称为路径的长度。

如果一个有向图从其中任何一个节点出发可以到达其他任何一个节点,则称这个有向图是强连通的。

如果从有向图的一个节点出发到返回这个节点的路径的长度都是kk的倍数(k>1,k∈Nk>1, k\in \mathcal{N}),则称这个节点是周期性节点。如果一个有向图不含周期性节点,则称这个有向图是非周期性图(Aperiodic Graph),否则为周期性图。

2、随机游走模型

给定一个含有nn个节点的有向图,在有向图上定义随机游走模型,也即一阶马尔科夫链。其中节点表示状态,有向边表示状态之间的转移。假设一个节点通过有向边到达其他节点的概率是相同的。

具体来说,转移矩阵M∈Rn×n\boldsymbol{M}\in \mathbb{R}^{n\times n},其中Mij\boldsymbol{M}_{ij}表示从节点jj到节点ii的概率。假设节点jj有kk条出边,对于jj连出的节点ii,有Mij=1/k\boldsymbol{M}_{ij}=1/k;对于其他节点,Mij=0\boldsymbol{M}_{ij}=0。

显然,转移矩阵满足以下性质:

  • mij⩾0m_{ij}\geqslant 0
  • ∑i=1nmij=1,∀j\sum_{i=1}^n m_{ij}=1, \forall j

转移矩阵就是一个随机矩阵。


随机游走者每经过一个单位时间转移一个状态,假设当前时刻在第jj个节点(状态为jj),下一个时刻转移到第ii个节点的概率为Mij\boldsymbol{M}_{ij}。显然,这一概率只依赖于当前状态,与过去的状态无关,具有马尔科夫性质,构成一个一阶马尔科夫链。

随机游走者在某个时刻tt访问各个节点的概率,可以用一个nn维列向量RtR_t表示,这也就是马尔科夫链在tt时刻的状态分布。则有:

Rt+1=MTRtR_{t+1}=\boldsymbol{M}^T R_t

3、PageRank 问题

将网页看作节点,网页之间的跳转看作边,构成一个有向图。浏览者随机浏览网页,构成一个随机游走模型。

PageRank 是一个函数,输入是网页,输出是一个实数值,表示这个网页的重要性。得到 PageRank 值的一种方法是:假定浏览者随机游走的情况下,考虑其某个时刻停留在某个页面的概率,这个概率值就作为这个页面的 PageRank 值。所以 PageRank 就是该随机游走模型的平稳分布,每个页面的 PageRank 值就是平稳概率。

3.1、基本定义

给定一个包含nn个节点v1,v2,…,vnv_1, v_2, \dots, v_n的强连通的,非周期性的有向图。在有向图上定义随机游走模型,即一阶马尔科夫链,其转移矩阵为M\boldsymbol{M}。这个马尔科夫链具有平稳分布R\boldsymbol{R},称其为这个有向图的 PageRank。R\boldsymbol{R}的各个分量称为各个节点的 PageRank 值,记为PR(vi)PR(v_i)。


考虑在时刻0,1,2,…,t,…0, 1, 2, \dots, t, \dots,访问各个节点的概率分布为R0,MR0,M2R0,…,MtR0,…\boldsymbol{R}_0, \boldsymbol{M}\boldsymbol{R}_0, \boldsymbol{M}^2\boldsymbol{R}_0, \dots, \boldsymbol{M}^t\boldsymbol{R}_0, \dots

由于不可约且非周期的有限状态马尔科夫链,由唯一平稳分布存在,并且当时间趋于无穷时的状态分布收敛于唯一的平稳分布。PageRank 问题的基本定义满足上述条件,所以极限:

lim⁡t→+∞MtR0=R\lim_{t\to +\infty} \boldsymbol{M}^t\boldsymbol{R}_0=\boldsymbol{R}

存在。极限值R\boldsymbol{R}就表示平稳分布,满足MR=R\boldsymbol{M}\boldsymbol{R}=\boldsymbol{R}。


记M(vi)M(v_i)为存在到viv_i的出边的节点的集合,L(vj)L(v_j)为节点vjv_j的出度。

PageRank 值满足以下性质:

  • PR(vi)⩾0PR(v_i)\geqslant 0
  • ∑i=1nPR(vi)=1\sum_{i=1}^n PR(v_i)=1
  • PR(vi)=∑vj∈M(vi)PR(vj)L(vj)PR(v_i)=\sum_{v_j\in M(v_i)}\frac{PR(v_j)}{L(v_j)}

PageRank 的定义是明确的,可以通过迭代的方式求出R\boldsymbol{R}。

3.2、一般定义

有时候有向图并不满足强连通和非周期性的条件,则此时不一定能够得到有意义的R\boldsymbol{R}。

这就引出了 PageRank 的一般定义,其基本思想是在基本定义的基础上加入平滑项。考虑另一个完全随机游走模型,其转移矩阵的元素全部为1/n1/n。两个转移矩阵的线性组合构成一个新的转移矩阵,可以证明,其对应的马尔科夫链存在平稳分布R\boldsymbol{R},满足

R=dMR+1−dn1\boldsymbol{R}=d \boldsymbol{M}\boldsymbol{R}+\frac{1-d}{n}\boldsymbol{1}

其中1\boldsymbol{1}为全 1 向量,d∈[0,1]d\in [0, 1]是阻尼因子(Damping Factor),取值由经验决定。当dd接近 1 时,表示随机游走主要依据原始转移矩阵进行;当dd接近 0 时,表示随机游走主要依据完全随机游走模型进行。


一般定义下的 PageRank 值满足以下性质:

PR(vi)=d(∑vj∈M(vi)PR(vj)L(vj))+1−dnPR(v_i)=d\left(\sum_{v_j\in M(v_i)}\frac{PR(v_j)}{L(v_j)}\right)+\frac{1-d}{n}

等式右边第二项极为平滑项,它保证了\PR(v_i)>0,且仍有∑i=1nPR(vi)=1\sum_{i=1}^n PR(v_i)=1。

3.3、PageRank 的计算

3.3.1、迭代法

有定义,直接按照以下迭代公式进行迭代即可,直到收敛或达到某种精度:

Rt+1=dMRt+1−dn1\boldsymbol{R}_{t+1}=d\boldsymbol{M}\boldsymbol{R}_t+\frac{1-d}{n}\boldsymbol{1}

3.3.2、幂法

考虑矩阵A∈Rn×n\boldsymbol{A}\in\mathbb{R}^{n\times n}。假设其有nn个特征值∣λ1∣>∣λ2∣⩾⋯⩾∣λn∣|\lambda_1|>|\lambda_2|\geqslant \dots \geqslant |\lambda_n|,对应的特征向量为u1,u2,…,un\boldsymbol{u}_1, \boldsymbol{u}_2, \dots, \boldsymbol{u}_n。称λ1\lambda_1为A\boldsymbol{A}的主特征值,u1\boldsymbol{u}_1为对应的主特征向量。幂法用于近似求解主特征向量。

nn个特征向量构成nn维空间的一组基。任取一个初始向量x0\boldsymbol{x_0},用这组基表示为x0=∑i=1nauix_0=\sum_{i=1}^n a \boldsymbol{u}_i。假设λ1\lambda_1是特征方程的单根,则有:

xk=Akx0=∑i=1naλikui=a1λ1k[u1+∑i=2n(λiλ1)kaiui]\boldsymbol{x_k}=\boldsymbol{A}^k x_0=\sum_{i=1}^n a \lambda_i^k \boldsymbol{u}_i=a_1 \lambda_1^k \left[\boldsymbol{u}_1+\sum_{i=2}^n \left(\frac{\lambda_i}{\lambda_1}\right)^k a_i \boldsymbol{u}_i\right]

所以有:

lim⁡k→+∞xk=a1λ1ku1\lim_{k\to +\infty}\boldsymbol{x_k}= a_1 \lambda_1^k \boldsymbol{u}_1

也即kk充分大时,xk\boldsymbol{x_k}与u1\boldsymbol{u}_1仅相差一个系数。


一般 PageRank 中,转移矩阵可写为:

R=dMR+1−dn1=(dM+1−dnE)R=AR\boldsymbol{R}=d \boldsymbol{M}\boldsymbol{R}+\frac{1-d}{n}\boldsymbol{1}=\left(d \boldsymbol{M}+\frac{1-d}{n}\boldsymbol{E}\right)\boldsymbol{R}=\boldsymbol{A}\boldsymbol{R}

其中E\boldsymbol{E}为全 1 矩阵(由于∑i=1nRi=1\sum_{i=1}^n R_i=1,所以有ER=1\boldsymbol{E}\boldsymbol{R}=\boldsymbol{1})

由 Perron-Frobenius 定理,R\boldsymbol{R}是A\boldsymbol{A}的主特征向量。所以,可以用幂法求出 PageRank。

为了防止计算过程中出现绝对值过大或者过小的情况,过程中对xk\boldsymbol{x}_k要进行规范化处理:xk←xk/∥xk∥\boldsymbol{x}_k\leftarrow\boldsymbol{x}_k/\|\boldsymbol{x}_k\|。

3.3.3、直接法

根据定义,可以直接写出R\boldsymbol{R}的表达式:

R=1−dn(I−dM)−11\boldsymbol{R}=\frac{1-d}{n}(\boldsymbol{I}-d\boldsymbol{M})^{-1}\boldsymbol{1}

求解一个逆矩阵即可。

4、最大流问题

4.1、定义

定义容量网络N=(V,E,c,s,t)N=(V, E, c, s, t),其中:

  • (V,E)(V, E)是有向连通图, 记n=∣V∣,m=∣E∣n=|V|, m=|E|
  • c:E→R+c:E\to \mathbb{R}^+是边的容量函数,表示边能够通过的最大流量
  • s,t∈Vs, t\in V是两个特殊的节点,分别是源点(发点)和汇点(收点),其余节点称为中间点。

流可以用一个函数f:E→R+f: E\to \mathbb{R}^+来表示,函数值表示边上的流量。若其满足下列条件(其中f((i,j))f((i, j))简记为f(i,j)f(i, j),函数cc同理):

  • 容量限制:f(i,j)⩽c(i,j),∀(i,j)∈Ef(i, j)\leqslant c(i, j), \forall (i, j)\in E,也即每条边上的流量不能超过其容量
  • 平衡条件:∑(j,i)∈Ef(j,i)−∑(i,j)∈Ef(i,j)=0,∀i∈V−{s,t}\sum_{(j, i) \in E} f(j, i)-\sum_{(i, j)\in E} f(i, j)=0, \forall i\in V-\{s, t\},也即中间点上流入量等于流出量

则称ff是网络NN上的一个可行流。


称源点ss的净流出量为ff的流量,记作v(f)v(f):

v(f)=∑(s,i)∈Ef(s,i)−∑(i,s)∈Ef(i,s)v(f)=\sum_{(s, i)\in E} f(s, i)-\sum_{(i, s)\in E} f(i, s)

流量最大的流称为最大流。显然有:

v(f)=∑(i,t)∈Ef(i,t)−∑(t,i)∈Ef(t,i)v(f)=\sum_{(i, t)\in E} f(i, t)-\sum_{(t, i)\in E} f(t, i)

也即源点的净流出量等于汇点的净流入量。


最大流问题的线性规划形式可以写为:

max⁡v(f)s.t.f(i,j)⩽c(i,j),∀(i,j)∈E∑(j,i)∈Ef(j,i)−∑(i,j)∈Ef(i,j)=0,∀i∈V−{s,t}v(f)=∑(s,i)∈Ef(s,i)−∑(i,s)∈Ef(i,s)f(i,j)⩾0,∀(i,j)∈Ev(f)⩾0\begin{align*} \max & v(f)\\ \text{s.t.} & f(i, j)\leqslant c(i, j), \forall (i, j)\in E\\ & \sum_{(j, i) \in E} f(j, i)-\sum_{(i, j)\in E} f(i, j)=0, \forall i\in V-\{s, t\}\\ & v(f)=\sum_{(s, i)\in E} f(s, i)-\sum_{(i, s)\in E} f(i, s) \\ & f(i, j)\geqslant 0, \forall (i, j)\in E\\ & v(f)\geqslant 0 \\ \end{align*}

能够求解线性规划的算法都能够求解最大流问题,但是最大流问题应用广泛,有很多专门的算法。

4.2、最小割集

设容量网络N=(V,E,c,s,t)N=(V, E, c, s, t),考虑将顶点集划分为两个集合,源点和汇点分别在两个集合中。也即A⊂VA\subset V且s∈A,t∈A‾s\in A, t\in \overline{A},则称:

(A,A‾)={(i,j)∈E∣i∈A,j∈A‾}(A, \overline{A})=\{(i, j)\in E|i\in A, j\in \overline{A}\}

为NN的一个割集。

扩展容量函数cc的定义,对整个割集求容量:

c(A,A‾)=∑(i,j)∈(A,A‾)c(i,j)c(A, \overline{A})=\sum_{(i, j)\in (A, \overline{A})} c(i, j)

容量最小的割集就称为最小割集。


引理 1:若ff是NN上的一个可行流,(A,A‾)(A, \overline{A})是一个割集,则有:

v(f)=∑(i,j)∈(A,A‾)f(i,j)−∑(j,i)∈(A,A‾)f(j,i)v(f)=\sum_{(i, j)\in (A, \overline{A})} f(i, j)-\sum_{(j, i)\in (A, \overline{A})} f(j, i)

也即整个网络上的流量,等于割集上的流量。

证明:从v(f)v(f)的定义出发,利用平衡条件,导出两者相等:

v(f)=∑(s,i)∈Ef(s,i)−∑(i,s)∈Ef(i,s)=∑(s,i)∈Ef(s,i)−∑(i,s)∈Ef(i,s)+∑i∈A−{s}(∑(j,i)∈Ef(j,i)−∑(i,j)∈Ef(i,j))=∑i∈A(∑(j,i)∈Ef(j,i)−∑(i,j)∈Ef(i,j))=∑i∈A∑(i,j)∈Ef(i,j)−∑i∈A∑(j,i)∈Ef(j,i)=(∑i∈A,j∈A∑(i,j)∈Ef(i,j)+∑i∈A,j∈A‾∑(i,j)∈Ef(i,j))−(∑i∈A,j∈A)∑(j,i)∈Ef(j,i)+∑i∈A,j∈A‾∑(j,i)∈Ef(j,i))=∑(i,j)∈(A,A‾)f(i,j)−∑(j,i)∈(A,A‾)f(j,i)\begin{align*} v(f)&=\sum_{(s, i)\in E} f(s, i)-\sum_{(i, s)\in E} f(i, s)\\ &=\sum_{(s, i)\in E} f(s, i)-\sum_{(i, s)\in E} f(i, s)+\sum_{i\in A-\{s\}}\left(\sum_{(j, i)\in E} f(j, i)-\sum_{(i, j)\in E} f(i, j)\right)\\ &=\sum_{i\in A}\left(\sum_{(j, i)\in E} f(j, i)-\sum_{(i, j)\in E} f(i, j)\right)\\ &=\sum_{i\in A}\sum_{(i, j)\in E} f(i, j)-\sum_{i\in A}\sum_{(j, i)\in E} f(j, i)\\ &=\left(\sum_{i\in A, j\in A}\sum_{(i, j)\in E} f(i, j)+\sum_{i\in A, j\in \overline{A}}\sum_{(i, j)\in E} f(i, j)\right)-\left(\sum_{i\in A, j\in A})\sum_{(j, i)\in E} f(j, i)+\sum_{i\in A, j\in \overline{A}}\sum_{(j, i)\in E} f(j, i)\right)\\ &=\sum_{(i, j)\in (A, \overline{A})} f(i, j)-\sum_{(j, i)\in (A, \overline{A})} f(j, i) \end{align*}


引理 2:若ff是NN上的一个可行流,(A,A‾)(A, \overline{A})是一个割集,则有:

v(f)⩽c(A,A‾)v(f)\leqslant c(A, \overline{A})

证明:

v(f)=∑(i,j)∈(A,A‾)f(i,j)−∑(j,i)∈(A,A‾)f(j,i)⩽∑(i,j)∈(A,A‾)f(i,j)⩽∑(i,j)∈(A,A‾)c(i,j)=c(A,A‾)\begin{align*} v(f)&=\sum_{(i, j)\in (A, \overline{A})} f(i, j)-\sum_{(j, i)\in (A, \overline{A})} f(j, i)\\ &\leqslant \sum_{(i, j)\in (A, \overline{A})} f(i, j)\\ & \leqslant \sum_{(i, j)\in (A, \overline{A})} c(i, j)=c(A, \overline{A}) \end{align*}


由引理 1、2,显然有引理 3:设ff是NN上的一个可行流,(A,A‾)(A, \overline{A})是一个割集。若v(f)=c(A,A‾)v(f)=c(A, \overline{A}),则ff是最大流,(A,A‾)(A, \overline{A})是最小割集。

4.3、最小割集的一个应用

经营投资策略问题。假设开发产品AiA_i需要先购进mim_i个工具Ti1,Ti2,…,TimiT_{i_1}, T_{i_2}, \dots, T_{i_{m_i}},工具TiT_i的价格为QiQ_i(每件工具只能用于一种产品的开发),产品AiA_i的收益为PiP_i。求选择哪些产品开发,使得利润最大。


考虑构建一个这样的容量网络:

capacity-net

如果产品AiA_i需要工具TjT_j,就拉一条从TjT_j到AiA_i的边,容量无限大的边。

以选择开发所有正净利润产品为例。将净利润为正的产品以及所需工具,放入汇点集合A‾;净利润为负的产品,\overline{A};净利润为负的产品,放入源点集合AA。分析割集的容量,割集中的边包括:

  • 从源点到正净利润产品所需的工具的边,这些边的容量之和代表成本
  • 从负净利润产品到汇点的边,这些边的容量之和代表负净利润产品的收益,我们不开发这些产品

用所有产品的收益之和减掉这个割集的容量,就是最终的利润。所以,求出最小割集的容量,就能得到最大利润。

4.4、增广链

设ff是NN上的一个可行流,定义以下概念:

  • 饱和边:流量等于容量的边
  • 非饱和边:流量小于容量的边
  • 零流边:流量为 0 的边
  • 非零流边:流量不为 0 的边

NN中从节点 ii 到jj的一条无重复边的路径(不考虑边的方向)称之为i-j 链。i-j 链的方向是从ii到jj,链中与链方向一致的边称为前向边,与链方向相反的边称为后向边。

i-j 增广链是前向边均为非饱和边,后向边均为非零流边的 i-j 链。


如果NN上存在一条关于可行流ff的 s-t 增广链,记Ef,EbE_f, E_b分别为增广链上的前向边和后向边,令:

δ=min⁡{min⁡e∈Ef(c(e)−f(e)),min⁡e∈Ebf(e)}\delta=\min\{\min_{e\in E_f} (c(e)-f(e)), \min_{e\in E_b} f(e)\}

可以构造一个新的流:

f′(e)={f(e)+δ,e∈Eff(e)−δ,e∈Ebf(e),otherwisef'(e)=\begin{cases} f(e)+\delta, & e\in E_f\\ f(e)-\delta, & e\in E_b\\ f(e), & otherwise \end{cases}

易有,f′f'也是一个可行流,且v(f′)=v(f)+δv(f')=v(f)+\delta,比原可行流的流量增加了δ\delta。


定理 1:可行流ff是最大流,当且仅当不存在关于ff的 s-t 增广链。

证明:由上面的叙述以及反证法,易证得必要性。下面证明充分性。

假设不存在关于ff的 s-t 增广链。考虑所有以ss为起点的增广链的终点jj的集合AA,由假设有t∉At\notin A。割集(A,A‾)(A, \overline{A})一定满足:

  • f(i,j)=c(i,j),∀(i,j)∈(A,A‾)f(i, j)=c(i, j), \forall (i, j)\in (A, \overline{A})
  • f(i,j)=0,∀(i,j)∈(A‾,A)f(i, j)=0, \forall (i, j)\in (\overline{A}, A)

对于其中第一点,若不然,存在(i,j)∈(A,A‾)(i, j)\in (A, \overline{A}),使得f(i,j)<c(i,j)f(i, j)<c(i, j),则 s-i 增广链延伸得到 s-j 链也是增广链,而 j∈A‾j\in \overline{A},矛盾。第二点,同理可证。

于是:

v(f)=∑(i,j)∈(A,A‾)f(i,j)−∑(j,i)∈(A,A‾)f(j,i)=∑(i,j)∈(A,A‾)c(i,j)=c(A,A‾)v(f)=\sum_{(i, j)\in (A, \overline{A})} f(i, j)-\sum_{(j, i)\in (A, \overline{A})} f(j, i)=\sum_{(i, j)\in (A, \overline{A})} c(i, j)=c(A, \overline{A})

由引理 3,可知ff是最大流。


定理 1 的证明中,还同时证明了 (A,A‾)(A, \overline{A})是一个最小割集,从而证明了以下定理:

定理 2(最大流-最小割定理):容量网络最大流的流量等于最小割集的容量。

5、Ford-Fulkerson 算法

简称 FF 算法,是解决最大流问题的一个经典算法。其基本思想是不断寻找增广链,修改链上的流量,直到不存在增广链为止,从而得到最大流。


从ss开始,逐个给顶点作标号,直到tt得到标号为止。某个顶点jj得到标号,表示已经找到从ss到jj的增广链,标号为(lj,δj)(l_j, \delta_j)。其中:

  • lj=+il_j=+i,表示增广链上jj的前一个顶点是ii,且(i,j)(i, j)是前向边;lj=−il_j=-i,表示增广链上jj的前一个顶点是ii,且(j,i)(j, i)是后向边
  • δj\delta_j是目前增广链上的δ\delta

将顶点分为三类:已标号已检查的、已标号未检查的TT、未标号的RR。FF 算法的完整流程叙述如下:

  1. 一开始,所有的点都是未标号的R=VR=V,边上的流量全为 0(零流)。
  2. 给ss标号(Δ,+∞)(\Delta, +\infin),T={s},R=V−{s}T=\{s\}, R=V-\{s\}
  3. 从TT中选取一个顶点ii
    1. 循环遍历RR中所有与ii邻接的所有顶点jj:
      • 如果 f(i,j)<c(i,j)f(i, j)<c(i, j),则给 jj 标号 (+i,δj)(+i, \delta_j)。其中 δj=min⁡{δi,c(i,j)−f(i,j)}\delta_j=\min\{\delta_i, c(i, j)-f(i, j)\}
      • 如果且 f(j,i)>0f(j, i)>0 ,则给 jj 标号 (−i,δj)(-i, \delta_j) 。其中 δj=min⁡{δi,f(j,i)}\delta_j=\min\{\delta_i, f(j, i)\}
    2. ii移出TT,jj从RR移入TT。
  4. 第 3 步中,若某轮循环中j=tj=t,说明找到了一条 s-t 增广链,该轮循环结束后,转入步骤 5;否则,回到步骤 2,直到TT为空。如果直到TT为空,都没有找到s−ts-t增广链,说明已经是最大流,结束算法。
  5. 从tt开始,沿着标号链回溯,修改流量。具体而言:

    {f(i,j)←f(i,j)+δj,j←i,lj=+if(j,i)←f(j,i)−δj,j←i,lj=−i\begin{cases} f(i, j)\leftarrow f(i, j)+\delta_j, j\leftarrow i, & l_j=+i\\ f(j, i)\leftarrow f(j, i)-\delta_j, j\leftarrow i, & l_j=-i \end{cases}

    直到回溯到ss,返回步骤 2。

一个应用 FF 算法求最大流的例子如下:

alt text

alt text


假设所有的容量都是正整数,记C=∑(i,j)∈Ec(i,j)C=\sum_{(i, j)\in E} c(i, j)。由于每次修改,流量至少增加 1,至多需要修改CC次。修改一次,至多需要O(m)O(m)的时间,所以 FF 算法的时间复杂度为O(mC)O(mC)。

计算机中,数的表示有一定的精度,所以算法总能在有限步内终止。理论上,容量为无理数时,δ\delta会越来越小,趋向于 0,算法不能在有限步终止。

6、Dinic 算法

FF 算法中,并没有明确给出标号过程的细节,找到的增广链以及找到不同增广链的顺序都是不固定的。而且 FF 算法每次只找一条增广链,就要重新标号,实际上也是一种浪费。

Dinic 算法对此进行了改进:

  • 每次求最短的 s-t 增广链
  • 充分利用一次标号的信息,每次找出尽可能多的增广链

6.1、辅助网络

定义一个关于容量网络NN和其上一个可行流ff的辅助网络 Nf=(V,Ef,ac,s,t)N_f=(V, E_f, ac, s, t),其中:

  • Ef+={(i,j)∣(i,j)∈E,f(i,j)<c(i,j)}E_f^+=\{(i, j)|(i, j)\in E, f(i, j)<c(i, j)\},Ef−={(j,i)∣(i,j)∈E,f(i,j)>0}E_f^-=\{(j, i)|(i, j)\in E, f(i, j)>0\}
  • Ef=Ef+∪Ef−E_f=E_f^+\cup E_f^-

acac是辅助容量函数,定义如下:

ac(i,j)={c(i,j)−f(i,j),(i,j)∈Ef+f(j,i),(i,j)∈Ef−ac(i, j)=\begin{cases} c(i, j)-f(i, j), & (i, j)\in E_f^+\\ f(j, i), & (i, j)\in E_f^-\\ \end{cases}

正向边的容量,是原网络上该边剩余的容量;反向边的容量,是原网络上该边当前的流量,也即能退回的容量。

显然NfN_f也是一个容量网络。这个容量网络也可以理解为一个“残差网络”,表示在原网络已经有的流 ff 之上,还能通过多少流量。

这样构建的一个辅助网络,是原网络的一个子网络。如果存在一条 s 到 t 的路径,且路径上每条边的辅助容量都大于 0,则这条路径就是原网络上的一条增广链。


引理 4:设NN的最大流量是v∗v^*,ff是NN上的一个可行流,则NfN_f的最大流量是v∗−v(f)v^*-v(f)(注意,辅助网络的定义与一个可行流有关,所以辅助网络的最大流当然与ff有关)。

证明:由于NN和NfN_f的点集相同。NN中的一个割集(A,A‾)(A, \overline{A}),NfN_f中也可以基于这两个集合定义出相应割集(A,A‾)f(A, \overline{A})_f,其容量:

ac(A,A‾)f=∑(i,j)∈Ef+[c(i,j)−f(i,j)]+∑(i,j)∈Ef−f(j,i)=∑(i,j)∈(A,A‾)[c(i,j)−f(i,j)]+∑(i,j)∈(A‾,A)f(j,i)=∑(i,j)∈(A,A‾)c(i,j)−[∑(i,j)∈(A,A‾)f(i,j)−∑(i,j)∈(A‾,A)f(j,i)]=c(A,A‾)−v(f)\begin{align*} ac(A, \overline{A})_f&=\sum_{(i, j)\in E_f^+}[c(i, j)-f(i, j)]+\sum_{(i, j)\in E_f^-}f(j, i)\\ &=\sum_{(i, j)\in (A, \overline{A})}[c(i, j)-f(i, j)]+\sum_{(i, j)\in (\overline{A}, A)}f(j, i)\\ &=\sum_{(i, j)\in (A, \overline{A})}c(i, j)-\left[\sum_{(i, j)\in (A, \overline{A})}f(i, j)-\sum_{(i, j)\in (\overline{A}, A)}f(j, i)\right]\\ &=c(A, \overline{A})-v(f) \end{align*}

由最大流最小割定理即得证。


设ff是NN上的一个可行流,gg是NfN_f上的一个可行流,补充定义g(i,j)=0,∀(i,j)∉Efg(i, j)=0, \forall (i, j)\notin E_f。再定义EE上的f′=f+gf'=f+g如下:

f′(i,j)=f(i,j)+g(i,j)−g(j,i),∀(i,j)∈Ef'(i, j)=f(i, j)+g(i, j)-g(j, i), \forall (i, j)\in E

则有引理 5:f′f'是NN上的一个可行流,且v(f′)=v(f)+v(g)v(f')=v(f)+v(g)。

证明:首先证明f′f'是一个可行流,要满足容量限制和平衡条件。

对于容量限制条件,由于0⩽g(i,j)⩽c(i,j)−f(i,j),0⩽g(j,i)⩽f(i,j)0\leqslant g(i, j)\leqslant c(i, j)-f(i, j), 0\leqslant g(j, i)\leqslant f(i, j),所以:

0⩽f′(i,j)=f(i,j)+g(i,j)−g(j,i)⩽c(i,j)0\leqslant f'(i, j)=f(i, j)+g(i, j)-g(j, i)\leqslant c(i, j)

得证。

对于平衡条件,只需证对于所有中间点ii有∑(j,i)∈Ef′(j,i)−∑(i,j)∈Ef′(i,j)=0\sum_{(j, i)\in E} f'(j, i)-\sum_{(i, j)\in E} f'(i, j)=0。∀i∈V−{s,t}\forall i\in V-\{s, t\},有:

∑(j,i)∈Ef′(j,i)=∑(j,i)∈Ef(j,i)+∑(j,i)∈E,(j,i)∈Efg(j,i)−∑(j,i)∈E,(i,j)∈Efg(i,j)∑(i,j)∈Ef′(i,j)=∑(i,j)∈Ef(i,j)+∑(i,j)∈E,(i,j)∈Efg(i,j)−∑(i,j)∈E,(j,i)∈Efg(j,i)\begin{align*} \sum_{(j, i)\in E} f'(j, i)&=\sum_{(j, i)\in E} f(j, i)+\sum_{(j, i)\in E, (j, i)\in E_f}g(j, i)-\sum_{(j, i)\in E, (i, j)\in E_f}g(i, j)\\ \sum_{(i, j)\in E} f'(i, j)&=\sum_{(i, j)\in E} f(i, j)+\sum_{(i, j)\in E, (i, j)\in E_f}g(i, j)-\sum_{(i, j)\in E, (j, i)\in E_f}g(j, i) \end{align*}

所以:

∑(j,i)∈Ef′(j,i)−∑(i,j)∈Ef′(i,j)=[∑(j,i)∈Ef(j,i)−∑(i,j)∈Ef(i,j)]+[∑(j,i)∈Efg(j,i)−∑(i,j)∈Efg(i,j)]=0+0=0\begin{align*} \sum_{(j, i)\in E} f'(j, i)-\sum_{(i, j)\in E} f'(i, j)&=\left[\sum_{(j, i)\in E} f(j, i)-\sum_{(i, j)\in E} f(i, j)\right]+\left[\sum_{(j, i)\in E_f}g(j, i)-\sum_{(i, j)\in E_f}g(i, j)\right]\\ &=0+0=0 \end{align*}

因此,f′f' 是 NN 上的一个可行流。

又因为

v(f′)=∑(s,j)∈E[f(s,j)+g(s,j)−g(j,s)]−∑(j,s)∈E[f(j,s)+g(j,s)−g(s,j)]=[∑(s,j)∈Ef(s,j)+∑(s,j)∈E,(s,j)∈Efg(s,j)−∑(s,j)∈E,(j,s)∈Efg(j,s)]−[∑(j,s)∈Ef(j,s)+∑(j,s)∈E,(j,s)∈Efg(j,s)−∑(j,s)∈E,(s,j)∈Efg(s,j)]=[∑(s,j)∈Ef(s,j)−∑(j,s)∈Ef(j,s)]+[∑(s,j)∈Efg(s,j)−∑(j,s)∈Efg(j,s)]=v(f)+v(g)\begin{align*} v(f')&=\sum_{(s, j)\in E}[f(s, j)+g(s, j)-g(j, s)]-\sum_{(j, s)\in E}[f(j, s)+g(j, s)-g(s, j)]\\ &=\left[ \sum_{(s, j)\in E}f(s, j)+\sum_{(s, j)\in E, (s, j)\in E_f}g(s, j)-\sum_{(s, j)\in E, (j, s)\in E_f}g(j, s) \right]-\left[ \sum_{(j, s)\in E}f(j, s)+\sum_{(j, s)\in E, (j, s)\in E_f}g(j, s)-\sum_{(j, s)\in E, (s, j)\in E_f}g(s, j) \right]\\ &=\left[ \sum_{(s, j)\in E}f(s, j)-\sum_{(j, s)\in E}f(j, s) \right] + \left[ \sum_{(s, j)\in E_f}g(s, j)-\sum_{(j, s)\in E_f}g(j, s) \right]\\ &=v(f)+v(g) \end{align*}

故引理得证。

这个引理提供了一个标准操作:当你在原始网络中有一个流 ff,你可以在它的残余网络 NfN_f 上再找一个流 gg(可以想象为 ff 留下的“可用空间”或“可撤销的流量”上走的新流),然后按照公式将 gg“叠加”回 ff,从而在 NN 上得到一个新的、流量更大的可行流。

对于 gg 在正向边的流量,也即 (i,j)∈E(i, j)\in E,直接叠加上去即可,所以 +g(i,j)+g(i, j);对于 gg 在反向边的流量,也即 (j,i)∈E(j, i)\in E,则相当于从原网络的 (i,j)(i, j) 边上撤销一部分流量,所以 −g(j,i)-g(j, i)。


有了引理 4 和 5,Dinic 算法的基本思路就很清晰了:通过找到 NfN_f 的最大流,来找到 NN 的最大流。

6.2、分层辅助网络

在NfN_f中,记顶点ii在以ss为起点的广度优先生成树中的层数为d(i)d(i)(ss在第 0 层)。定义NfN_f的分层辅助网络ANf=(Vf,AEf,ac,s,t)AN_{f}=(V_f, AE_f, ac, s, t),其中:

  • Vf(k)={i∈Vf∣d(i)=k},0⩽k⩽d(t)V_{f}^{(k)}=\{i\in V_f|d(i)=k\}, 0\leqslant k\leqslant d(t)
  • Vf=⋃k=0d(t)Vf(k)V_f=\bigcup_{k=0}^{d(t)} V_f^{(k)}

分层辅助网络的边集定义为:

AEf=⋃k=0d(t)−1{(i,j)∈Ef∣i∈Vf(k),j∈Vf(k+1)}AE_f=\bigcup_{k=0}^{d(t)-1}\{(i, j)\in E_f|i\in V_f^{(k)}, j\in V_f^{(k+1)}\}

以下是一个带可行流的原始流量网络以及相应的辅助网络和分层辅助网络的示例,可以帮助理解三个网络的定义以及之间的关系:

multi-net

可以直观地理解为,多层辅助网络是按照广度优先搜索顺序进行搜索时,看到的网络。


通过建立分层辅助网络,能保证每次都能找到最短的增广链。而且由于分层辅助网络是一个有向无环图,避免了走到环上浪费时间。

6.3、Dinic 最大流算法

定义分层辅助网络中,每一个顶点ii的流通量th(i)th(i):

th(i)=min⁡{∑(j,i)∈AEfac(j,i),∑(i,j)∈AEfac(i,j)}th(i)=\min\{\sum_{(j, i)\in AE_f} ac(j, i), \sum_{(i, j)\in AE_f} ac(i, j)\}

规定 ac(i,s)=+∞,ac(t,j)=+∞,∀i∈AEfac(i, s)=+\infty, ac(t, j)=+\infty, \forall i\in AE_f

则 Dinic 最大流算法的完整叙述如下:

  1. 对于原始网络NN以及零流ff,构建AN(f)AN(f)
  2. 如果AN(f)AN(f)不存在从ss到tt的路径,返回当前ff,这就是NN的最大流,结束算法;否则进入步骤 3
  3. 初始化流gg为零流。
  4. 查看ANfAN_f中,流通量为 0 的顶点。如果其中包含ss或者tt,直接令f←f+gf\leftarrow f+g,返回步骤 2;否则,删去这些顶点以及与其相连的边。
  5. 找到流通量最小的顶点kk,从kk开始,将th(k)th(k)个单位的流往两边推送分别到达ss和tt(由流通量的定义,如果流要分流给多条边,一定是每条边都恰好占满),同时更新gg。然后删掉kk以及与其相连的边。(如果从某个顶点向外推送流的时候,有多条边可去,优先给容量较小的边)。
  6. 对于ANfAN_f中的每一条边,如果g(i,j)=ac(i,j)g(i, j)=ac(i, j),则删去这条边;否则ac(i,j)←ac(i,j)−g(i,j)ac(i, j)\leftarrow ac(i, j)-g(i, j),回到步骤 4。

对于上面那个示例中的分层辅助网络,Dinic 算法的执行过程如下:

  1. 首先可以依次删掉节点 9、7、5,因为它们的流通量会(变为) 0。然后找到流通量最小的顶点 1,其流通量为 3,向两边推送。

    alt text

  2. 推送完毕后,删掉节点 1 和与其相连的边,再更新辅助网络的边容量。

    alt text

  3. 接下来,继续找流通量最小的顶点。顶点 4 的流通量最小,为 2,重复上面的操作。会发现已经不存在从 ss 到 tt 的路径了,这一轮结束。此时就找到了极大流,更新到 ff 中。

    alt text

  4. 根据新的 ff,重新构建分层辅助网络,继续上述过程,直到整个算法结束。

6.4、Dinic 算法的细节

引理 6:在 Dinic 算法的某一轮结束时,分层辅助网络 ANfAN_f 中不存在从 ss 到tt的路径,此时 gg 是 ANfAN_f 的极大流。

证明:在 ANfAN_f 中, 只有当流通量变成 0 后, 才删去这个顶点及其
关联的边。这样的顶点和边不可能出现在关于 gg 的 s-t 前向增
广链上。当每个阶段结束时, ANfAN_f 剩余的图中已不存在从 s 到
t 的路径, 从而得到的 gg 是 ANfAN_f 的极大流。

这个引理保证了 Dinic 算法每一轮结束时,找到了分层辅助网络的极大流,从而保证了算法的正确性。


引理 7:在每个阶段, ANf+gAN_{f+g} 中 s-t 距离大于前一个阶段 ANfAN_f
中 s-t 距离(约定不存在 s-t 最短路径时, 视其距离为 +∞+\infty)

证明:记 f′=f+gf'=f+g。证明的核心是通过分析在加上 gg 前后边的变化情况,来推断最短路径的变化。

首先考虑原来在 ANfAN_f 中的边,是否会存在于 Nf+gN_{f+g} 中。

对于 ∀(i,j)∈AEf\forall (i, j)\in AE_f(一定有 d(i)<d(j)d(i)<d(j)),若同时 (i,j)∈E(i, j)\in E(是原网络中的边,也即是辅助网络中的前向边 Ef+E_f^+),则 ac(i,j)=c(i,j)−f(i,j)ac(i, j)=c(i, j)-f(i, j),又因为 f′(i,j)=f(i,j)+g(i,j)f'(i, j)=f(i, j)+g(i, j),所以:

  • 当 g(i,j)<ac(i,j)g(i, j)<ac(i, j) 时,有 f′(i,j)=c(i,j)−[ac(i,j)−g(i,j)]<c(i,j)f'(i, j)=c(i, j)-[ac(i, j)-g(i, j)]<c(i, j)。(i,j)(i, j) 未饱和,会出现在下一轮中,也即 (i,j)∈Ef+g(i, j)\in E_{f+g}

  • 当 g(i,j)=ac(i,j)g(i, j)=ac(i, j) 时,有 f′(i,j)=c(i,j)f'(i, j)=c(i, j),(i,j)(i, j) 达到饱和,不会再出现在下一轮中,也即 (i,j)∉Ef+g(i, j)\notin E_{f+g}

  • 此外,当 f(i,j)=0f(i, j)=0 且 g(i,j)>0g(i, j)>0,则 Nf+gN_{f+g} 会比 NfN_f 多出一条边 (j,i)(j, i)。因为 d(j)>d(i)d(j)>d(i),所以这条边是由高层到低层的边,不会出现在 ANf+gAN_{f+g} 中。

对于 ∀(j,i)∈Ef\forall (j, i)\in E_f(一定有 d(j)<d(i)d(j)<d(i)),若同时 (i,j)∈E(i, j)\in E(是原网络中的边,也即是辅助网络中的后向边 Ef−E_f^-),则 ac(j,i)=f(i,j)ac(j, i)=f(i, j),又因为 f′(i,j)=f(i,j)−g(j,i)f'(i, j)=f(i, j)-g(j, i),所以

  • 当 g(j,i)<ac(j,i)g(j, i)<ac(j, i) 时,有 f′(i,j)=ac(i,j)−g(j,i)>0f'(i, j)=ac(i, j)-g(j, i)>0。(i,j)(i, j) 为非零流边,(j,i)(j, i)会出现在下一轮中,也即 (j,i)∈Ef+g(j, i)\in E_{f+g}

  • 当 g(j,i)=ac(j,i)g(j, i)=ac(j, i) 时,有 f′(i,j)=0f'(i, j)=0,(i,j)(i, j) 变为零流边,(j,i)(j, i)不会再出现在下一轮中,也即 (j,i)∉Ef+g(j, i)\notin E_{f+g}

  • 此外,当 f(i,j)=c(i,j)f(i, j)=c(i, j) 但 f′(i,j)<c(i,j)f'(i, j)<c(i, j),则 Nf+gN_{f+g} 会比 NfN_f 多出一条边 (i,j)(i, j)。因为 d(i)>d(j)d(i)>d(j),所以这条边是由高层到低层的边,不会出现在 ANf+gAN_{f+g} 中。

综上所述,ANfAN_f 中的边如果也是 Nf+gN_{f+g} 的边,当且仅当它是关于 gg 的非饱和边。

因为加上 gg 产生的,在 Nf+gN_{f+g} 中而不在 ANfAN_f 中的新边,都是由高层到低层的边,不会出现在 ANf+gAN_{f+g} 中。其余的在 NfN_f 中而不在 ANfAN_f 中的边,都与 gg 无关,会直接保留到 Nf+gN_{f+g} 中。而且因为不在 ANfAN_f 中,这些边一定不是从低层指向高层的边。总之,在 Nf+gN_{f+g} 而不在 ANfAN_f 中的边,都不可能是低层指向高层的边。

由 ANf+gAN_{f+g} 本身构造时的性质,ANf+gAN_{f+g} 中从 s 到 t 的最短路径就是 Nf+gN_{f+g} 中从 s 到 t 的最短路径。

假设最短路径上的边全部由原本 ANfAN_f 中的边构成,由前所述,这些边一定都是关于 gg 的非饱和边。这意味着这条路径也是 ANfAN_f 中的一条路径,且路径上的边还有剩余容量,这条路径就是 ANfAN_f 中的一条 s-t 增广链,这与 gg 是 ANfAN_f 的极大流矛盾。

所以最短路径上至少有一条边不是原本 ANfAN_f 中的边。由前所述,这条边一定不是由低层指向高层的边,所以一定“绕了路”,从而使得这条路径的长度增加了至少 1。

因此,ANf+gAN_{f+g} 中 s-t 距离大于前一个阶段 ANfAN_f 中 s-t 距离,证毕。


定理 3:Dinic 算法最终得到的 ff 是最大流,且算法在 O(n3)O(n^3) 步内终止。

证明:当算法终止时,ANfAN_f 中不存在从 ss 到 tt 的路径,由 ANfAN_f 的定义,NfN_f 中也不存在从 ss 到 tt 的路径。所以,NfN_f 的最大流是 0,最大流量为 0。由引理 4,可知 v(f)=v∗v(f)=v^*,所以 ff 是最大流。

由引理 7,每个阶段结束时,ANfAN_f 中 s-t 距离至少增加 1。由于 s-t 距离最多为 n−1n-1,所以 Dinic 算法最多进行 n−1n-1 个阶段。

在每个阶段:

  • 构造 ANfAN_f —— O(m)O(m)
  • 计算顶点的流通量 —— O(m)O(m)
  • 找到流通量最小的顶点 —— O(n)O(n),至多需要找 O(n)O(n) 次 —— O(n2)O(n^2)
  • 从流通量最小的顶点向外推送流。
    • 第一种情况是,推流之后变成饱和边,删去改变,最多有 O(m)O(m) 次。
    • 第二种情况是,推流之后变成非饱和边,由于优先给容量较小的边,每次从某个顶点向外推送流时,至多产生一条非饱和边,故每次推流,最多有 O(n)O(n) 次。至多进行 nn 轮,共 O(n2)O(n^2) 次。

综上所述,每个阶段的时间复杂度为 O(n2+m)O(n^2 + m),整个算法的时间复杂度为 O(n3+nm)O(n^3 + nm),即 O(n3)O(n^3)。

6.5、简单容量网络上的 Dinic 算法

边的容量均为 1,中间点或入度为一或出度为一的容量网络,称为简单容量网络。

引理 8:在简单容量网络中,s 到 t 的距离 l⩽n/v∗+1l\leqslant n/v^* +1

证明:记到 s 的举例为 ii 的顶点集为 Vi,1⩽i⩽l−1V_i, 1\leqslant i \leqslant l-1。因为任何最大流一定要经过这些顶点,而任何中间顶点的入度或出度为 1,所以有 v∗⩽∣Vi∣v^* \leqslant |V_i|,所以有:

(l−1)v∗⩽∑i=1l−1∣Vi∣⩽n(l-1)v^* \leqslant \sum_{i=1}^{l-1} |V_i| \leqslant n

所以有 l⩽n/v∗+1l\leqslant n/v^* +1。

定理 4:简单容量网络上,Dinic 算法在O(n1/2m)O(n^{1/2}m)步内终止。

7、最大流的应用

7.1、不相交路径问题

Edge Disjoint Paths,不相交路径。有向图上,对于给定的两个顶点,如果这两个顶点之间的两条路径没有公共边,就称这两条路径是不相交的。对于给定的两个顶点,求可以同时存在的不相交路径的最大条数。


考虑将给定的两个顶点看作源点和汇点,所有边的容量都设为 1(保证只能有 1 个单位的流流过,也即只能走一次。在这样的容量网络上求最大流,最大流的流量即为答案。


对于由多个起点和多个终点的不相交路径问题。考虑额外引入一个超级源点和一个超级汇点,超级源点到每个起点连一条容量无限的边,每个终点到超级汇点也连一条容量无限的的边。在这样的容量网络上求最大流,就可以得出结果了。

7.2、顶点存在容量限制的最大流

假设顶点也有容量,表示最多有多少流量能够流入该节点。

假设顶点vv有容量限制c(v)c(v),考虑创建一个新顶点v′v',从vv到v′v'连一条容量为c(v)c(v)的边,然后将所有原来vv的出边改为从v′v'出发。这样,就加个顶点容量转化为了边的容量,变成了普通的最大流问题。

在新网络上求解出最大流之后,将 vv 和 v′v' 合并回去就得到了原网络上的解。

7.3、独立路径问题

起点和终点相同的两条路径,如果其没有其他的公共顶点,就称这两条路径是独立的。对于给定的起点和终点,求可以同时存在的独立路径的最大条数。


给所有的中间点设置 1 的顶点容量限制。由 7.2 的方法,转化为普通最大流问题求解即可。

8、Floyd-Warshall 算法

也常简称为 Floyd 算法,这是一种在带负权边的图上求最短路径和检测负回路的算法

带权有向图D=(V,E,w)D=(V, E, w),其中w:E→Rw: E\rightarrow R是权函数,函数值代表边的权值。DD中,边的权值和为负的贿赂称之为负回路。

DD中任意两点之间要么有最短路径,要么不存在路径当且仅当DD中不存在负回路。


记d(k)(i,j)d^{(k)}(i, j)为从ii到jj,只经过{1,2,…,k}\{1, 2, \dots, k\}的顶点的最短路径长度。则有递推关系:

d(k)(i,j)={w(i,j),k=0min⁡{d(k−1)(i,j),d(k−1)(i,k)+d(k−1)(k,j)},1⩽k⩽n,1⩽i,j⩽n,i≠k,j≠kd^{(k)}(i, j)=\begin{cases} w(i, j), & k=0\\ \min\{d^{(k-1)}(i, j), d^{(k-1)}(i, k)+d^{(k-1)}(k, j)\}, & 1\leqslant k\leqslant n, 1\leqslant i, j\leqslant n, i\neq k, j\neq k \end{cases}

规定w(i,i)=0w(i, i)=0,且对(i,j)∉E(i, j)\notin E,w(i,j)=+∞w(i, j)=+\infin。

若DD中存在负回路,假设负回路经过ii,回路中除ii外顶点的最大号码为kk,则一定有d(k)(i,i)<0d^{(k)}(i, i)<0。


Floyd-Warshall 算法的时间复杂度为O(n3)O(n^3)。

9、最小费用流

9.1、定义

在容量网络N=(V,E,c,s,t)N=(V, E, c, s, t)的基础上,添加费用函数w:E→Rw: E\rightarrow R,表示每条边的单位流量费用,得到容量-费用网络N=(V,E,c,w,s,t)N=(V, E, c, w, s, t)。

设ff是NN上的一个可行流,称w(f)=∑(i,j)∈Ew(i,j)f(i,j)w(f)=\sum_{(i, j)\in E} w(i, j)f(i, j)为ff的费用。所有流量为vv的可行流中,费用最小的称为流量vv的最小费用流。

最小费用流问题就是,在容量-费用网络NN上,求流量为v0v_0的最小费用流。


容量-费用网络上也可以定义关于可行流ff的辅助网络Nf=(V,Ef,ac,aw,s,t)N_f=(V, E_f, ac, aw, s, t)。其中EfE_f和acac的定义与之前相同,辅助费用的定义如下:

aw(i,j)={w(i,j),(i,j)∈Ef+−w(j,i),(i,j)∈Ef−aw(i, j)=\begin{cases} w(i, j), & (i, j)\in E_f^+\\ -w(j, i), & (i, j)\in E_f^- \end{cases}

也有类似的引理 9:设ff是容量-费用网络NN上的可行流,gg是辅助网络NfN_f上的可行流,则f′=f+gf'=f+g也是NN上的一个可行流,且w(f′)=w(f)+aw(g)w(f')=w(f)+aw(g)。

9.2、圈流

设CC是容量-费用网络NN中一条边不重复的回路,ECE_C是CC的边集,CC上的圈流hCh^C定义为:

hC(i,j)={δ,(i,j)∈EC0,otherwiseh^C(i, j)=\begin{cases} \delta, & (i, j)\in E_C\\ 0, & otherwise \end{cases}

其中δ=min⁡(i,j)∈ECc(i,j)>0\delta = \min_{(i, j)\in E_C} c(i, j)>0称之为hCh^C的环流量。

显然,圈流是可行流,且v(hC)=0,w(hC)=δ∑(i,j)∈ECw(i,j)v(h^C)=0, w(h^C)=\delta\sum_{(i, j)\in E_C} w(i, j)。简记w(C)=∑(i,j)∈ECw(i,j)w(C)=\sum_{(i, j)\in E_C} w(i, j),aw(C)aw(C)同理。


设ff是NN上的一个可行流,由引理,f′=f+hCf'=f+h^C也是一个可行流,且v(f′)=v(f),w(f′)=w(f)+δ⋅aw(C)v(f')=v(f), w(f')=w(f)+\delta\cdot aw(C)。

可以这样求最小费用流:首先求一个流量为v0v_0的可行流。如果NfN_f中存在aw(C)<0aw(C)<0的负回路CC,则可以令f′←f+hCf'\leftarrow f+h^C(流量不变,费用减小),重复这个过程,直到NfN_f中不存在负回路为止。


例如,对下面这个容量-费用网络,构建出辅助网络之后,可以找到一个圈流。可以确定 δ=3\delta=3,然后更新流量:

alt text

9.3、最小费用流的负回路算法

设 f1,f2f_1, f_2 是两个可行流,定义在 Nf2N_{f_2} 上的流 g=f1−f2g=f_1-f_2 如下,∀(i,j)∈E\forall (i, j)\in E:

  • 若 f1(i,j)<f2(i,j)f_1(i, j) < f_2(i, j),则 g(i,j)=f2(i,j)−f1(i,j)g(i, j) = f_2(i, j) - f_1(i, j)
  • 若 f1(i,j)>f2(i,j)f_1(i, j) > f_2(i, j),则 g(i,j)=f1(i,j)−f2(i,j)g(i, j) = f_1(i, j) - f_2(i, j)
  • Nf2N_{f_2} 上其他边的 gg 值均为 0

有引理 10:设 f1,f2f_1, f_2 是两个可行流,g=f1−f2g=f_1-f_2 则 gg 是 Nf2N_{f_2} 上的一个可行流,且 f1=f2+gf_1=f_2+g。

可类似引理 5 进行证明。


引理 11:若可行流 ff 不是零流且 v(f)=0v(f)=0,则 ff 等于若干个圈流之和。

证明:记 E1={(i,j)∈E∣f(i,j)>0}E_1=\{(i, j)\in E | f(i, j) > 0\},D1D_1 是由 E1E_1 中的边以及相关联的顶点构成的子图。由于 ff 不是零流,所以 E1≠∅E_1 \neq \emptyset。

由于 ff 是可行流,且 D1D_1 是原网络的子图,根据平衡条件,对 D1D_1 中的任意顶点 ii,其在原网络中,流入都等于流出。假设在 D1D_1 中,ii 的出度为 0,由 D1D_1 的定义,这意味这原网络中,不存在 ii 出发的某条边 (i,k)(i, k) 使得 f(i,k)>0f(i, k)>0,则 ii 的流出为 0。由于 ff 是可行流,ii 的流入也为 0,则 E1E_1 中不存在与 ii 关联的边,这与 ii 是 D1D_1 的顶点矛盾。故 D1D_1 中每个顶点的出度至少为 1。

因为每个顶点的出度至少为 1,从任何一个顶点出发,沿着出边走下去,由于顶点数有限,一定会回到某个已经经过的顶点(鸽巢原理),从而形成一个回路,记这个回路为 CC。且一定可以找出一个不含重复边的回路(到第一个重复到达的顶点就停止)。令 δ1=min⁡{f(i,j)∣(i,j)∈EC}>0\delta_1=\min\{f(i, j)|(i, j)\in E_C\}>0,则 hC1h^{C_1} 是 C1C_1 上环流量为 δ1\delta_1 的圈流。

令 f1=f−hC1f_1 = f-h^{C_1},则易证得 f1f_1 也是一个可行流,且 v(f1)=0v(f_1)=0,而且 f1f_1 至少比 ff 多一条零流边。由于 ff 只有有限条零流边,重复上述过程,可以剥离出若干个圈流,直到剩下的流为零流为止。证毕。


定理 5:设 ff 是 NN 上流量为 v0v_0 的可行流,则 ff 是最小费用流当且仅当 NfN_f 中不存在以辅助费用 awaw 为权的负回路。

证明:之前已经说明了必要性,现证充分性。

假设 NfN_f 中不存在关于 awaw 的负回路。设 f′f' 是 NN 上任一流量为 v0v_0 的可行流,令 g=f′−fg=f'-f。由引理 10,gg 是 NfN_f 上的一个可行流,且 f′=f+gf'=f+g。又由引理 5,因为 v(f′)=v(f)v(f')=v(f),所以 v(g)=0v(g)=0。

不妨设 f′≠ff'\neq f,则 gg 不是零流,由引理 11,gg 可以表示为若干个圈流之和:

g=hC1+hC2+⋯+hCkg=h^{C_1}+h^{C_2}+\dots+h^{C_k}

所以 gg 的辅助费用:

aw(g)=δ1aw(C1)+δ2aw(C2)+⋯+δkaw(Ck)aw(g)=\delta_1 aw(C_1)+\delta_2 aw(C_2)+\dots+\delta_k aw(C_k)

其中由 δi>0\delta_i>0。由于 NfN_f 中不存在负回路,所以 ∀i,aw(Ci)⩾0\forall i, aw(C_i)\geqslant 0,从而 aw(g)⩾0aw(g)\geqslant 0。又由引理 9,w(f′)=w(f)+aw(g)⩾w(f)w(f')=w(f)+aw(g)\geqslant w(f),即 ff 是最小费用流,证毕。


最小费用流的负回路的完整算法叙述如下:

  1. 调用最大流算法。若求出的最大流流量小于v0v_0,无解,结束算法
  2. 构造NfN_f
  3. 用 Floyd 算法检测NfN_f中是否存在负回路,边的权函数是awaw。若不存在负回路,返回当前ff,结束算法
  4. 令 δ=min⁡{ac(i,j)∣(i,j)∈EC}\delta=\min\{ac(i, j)|(i, j)\in E_C\},然后更新流。具体而言:

    {f(i,j)←f(i,j)+δ,(i,j)∈ECf(j,i)←f(j,i)−δ,(j,i)∈EC\begin{cases} f(i, j)\leftarrow f(i, j)+\delta, & (i, j)\in E_C\\ f(j, i)\leftarrow f(j, i)-\delta, & (j, i)\in E_C \end{cases}

    然后回到步骤 2

9.4、最小费用流的最短路径算法

另一种求最小费用流的思路是,从一个某个费用最小的初始流ff开始(v(f)<v0v(f)<v_0),然后寻找一条费用最少的 s-t 增广链,修改PP上的流量,得到新的可行流f′f'。重复这个过程,直到流量达到v0v_0。关键是操作要保证 f′f' 是费用最小的。


引理 12:设有向图 D=(V,E)D=(V, E) 没有孤立点,顶点 ss 的出度比入度大 1,顶点 tt 的入度比出度大 1,其余顶点的入度等于出度。则 DD 可以表示为一条 s-t 路径与若干条回路的并。

证明:从 ss 开始,沿一条边到另一个顶点。由于除 s,ts, t 外的顶点的入度等于出度,到一个顶点后总可从一条没走过的边到另一个顶点, 直到 tt 为止. 这样得到一条 s-t 路径 P。

删去 P 涉及的边,此时可能出现一些孤立点,再把这些孤立点删去。则剩下的顶点,入度都等于出度且不为 0,它们由若干条回路组成。证毕。


定理 6:设 ff 是 NN 上的流量为 v0v_0 的最小费用流,PP 是 NfN_f 上的一条从 ss 到 tt 的最短路径(权值函数为 awaw),gg 是 PP 上流量为 θ\theta 的可行流,则 f′=f+gf'=f+g 是 NN 上流量为 v0+θv_0+\theta 的最小费用流。

证明:假设 f′f' 不是最小费用流。由定理 5,Nf′N_{f'} 中存在关于 awaw 的负回路 CC。

由于 f′f' 和 ff 仅在 PP 对应的增广链上不同, NfN_f 和 Nf′N_{f'} 边集之间的不同,一定与 PP 上的边有关。因此,若 (i,j)∈Nf′−Nf(i, j)\in N_{f'}-N_f,是新出现的边,则一定是在反向边推送流量导致的,也即 (j,i)∈EP(j, i)\in E_P。

因为 ff 是最小费用流,所以 NfN_f 中不存在关于 awaw 的负回路,所以 CC 中肯定包含不属于 NfN_f 的边。记这些边为 (i1,j1),(i2,j2),…,(ir,jr)(i_1, j_1), (i_2, j_2), \dots, (i_r, j_r)。则由前说明 (j1,i1),(j2,i2),…,(jr,ir)(j_1, i_1), (j_2, i_2), \dots, (j_r, i_r) 都在 PP 上。记这 2r2r 条边构成的集合为 HH。

从 PP 和 CC 一起构成的子图中删去 HH 以及因此产生的孤立点(PP 和 CC 中相同的边同时保留,视为两条一模一样的“平行边”),剩下的图记为 DD。观察 DD 可以发现,其满足引理 12 的条件,所以其可以视为一条从 ss 到 tt 的路径 P′P' 与若干条回路 C1,C2,…,ClC_1, C_2, \dots, C_l 的并。

于是有:

aw(D)=aw(P′)+aw(C1)+aw(C2)+⋯+aw(Cl)aw(D)=aw(P')+aw(C_1)+aw(C_2)+\dots+aw(C_l)

因为 aw(H)=0aw(H)=0 ,所以有:

aw(P)+aw(C)=aw(D)+aw(H)=aw(P′)+aw(C1)+aw(C2)+⋯+aw(Cl)aw(P)+aw(C)=aw(D)+aw(H)=aw(P')+aw(C_1)+aw(C_2)+\dots+aw(C_l)

所以有:

aw(P′)=aw(P)+aw(C)−[aw(C1)+aw(C2)+⋯+aw(Cl)]aw(P')=aw(P)+aw(C)-[aw(C_1)+aw(C_2)+\dots+aw(C_l)]

由于 C1,C2,…,ClC_1, C_2, \dots, C_l 都是 NfN_f 上的回路,而 ff 是最小费用流,所以 aw(Ci)⩾0,1⩽i⩽law(C_i)\geqslant 0, 1\leqslant i \leqslant l。又因为 aw(C)<0aw(C)<0,所以有 aw(P′)<aw(P)aw(P')<aw(P),这与 PP 是 NfN_f 上的权为 awaw 的最短路径矛盾,证毕。


最小费用流的最短路径算法的完整叙述如下:

  1. 初始化零流ff,当前流量v=0v=0
  2. 构造NfN_f
  3. 调用 Floyd 算法,计算NfN_f中ss到tt的最短路径,权值函数为awaw。如果不存在这样的路径,无解,结束算法;否则,记这条最短路径为PP
  4. 令 θ=min⁡{v,min⁡{ac(i,j)∣(i,j)∈EP}}\theta=\min\{v, \min\{ac(i, j)|(i, j)\in E_P\}\},然后更新流。具体而言,对EPE_P中的每一条边(i,j)(i, j):

    {f(i,j)←f(i,j)+θ,(i,j)∈Ef(j,i)←f(j,i)−θ,(j,i)∈E\begin{cases} f(i, j)\leftarrow f(i, j)+\theta, & (i, j)\in E\\ f(j, i)\leftarrow f(j, i)-\theta, & (j, i)\in E \end{cases}

  5. 令v←v+θv\leftarrow v+\theta,如果v⩾v0v\geqslant v_0,结束算法;否则,回到步骤 2

10、运输问题

10.1、问题定义与建模

有 mm 个产地 A1,A2,…,AmA_1, A_2, \dots, A_m,每个产地的产量为 aia_i;有 nn 个销地 B1,B2,…,BnB_1, B_2, \dots, B_n,每个销地的销量为 bjb_j。

已知从产地 AiA_i 到销地 BjB_j 运输 1 单位货物的费用为 wi,jw_{i,j}。要求设计一个运输方案,使得满足各产地的供货量和各销地的需货量的同时,总运输费用最小。


事实上,以上问题描述蕴含了一个产销平衡的假设,也即:

∑i=1mai=∑j=1nbj\sum_{i=1}^m a_i = \sum_{j=1}^n b_j

若产大于销,可以引入虚拟销地,需求量 bn+1=∑i=1mai−∑j=1nbjb_{n+1} = \sum_{i=1}^m a_i - \sum_{j=1}^n b_j,且运输费用 w(i,n+1)=0w(i, n+1) = 0。最终供给 bn+1b_{n+1} 的货物不需要实际运输出去,是剩余的产品量。

若销大于产,可以引入虚拟产地,产量 am+1=∑j=1nbj−∑i=1maia_{m+1} = \sum_{j=1}^n b_j - \sum_{i=1}^m a_i,且运输费用 w(m+1,j)=0w(m+1, j) = 0。最终 am+1a_{m+1} 供出的货物不需要实际运输出去,是缺少的产品量。

因此我们总是可以假设 ∑i=1mai=∑j=1nbj\sum_{i=1}^m a_i = \sum_{j=1}^n b_j。


运输问题可以看作一个最小费用流问题。构造一个容量-费用网络 N=(V,E,c,w,s,t)N=(V, E, c, w, s, t),其中:

  • V={s,t}∪{A1,A2,…,Am}∪{B1,B2,…,Bn}V=\{s, t\}\cup \{A_1, A_2, \dots, A_m\} \cup \{B_1, B_2, \dots, B_n\}
  • E={(s,Ai)∣1⩽i⩽m}∪{(Ai,Bj)∣1⩽i⩽m,1⩽j⩽n}∪{(Bj,t)∣1⩽j⩽n}E=\{(s, A_i)|1\leqslant i \leqslant m\} \cup \{(A_i, B_j)|1\leqslant i \leqslant m, 1\leqslant j \leqslant n\} \cup \{(B_j, t)|1\leqslant j \leqslant n\}
  • c(s,Ai)=ai,c(Ai,Bj)=+∞,c(Bj,t)=bjc(s, A_i)=a_i, c(A_i, B_j)=+\infin, c(B_j, t)=b_j
  • w(s,Ai)=0,w(Ai,Bj)=wi,j,w(Bj,t)=0w(s, A_i)=0, w(A_i, B_j)=w_{i, j}, w(B_j, t)=0

流量为 v0=∑i=1maiv_0=\sum_{i=1}^m a_i 的可行流就是运输问题的一个合法方案,而流量为 v0v_0 的最小费用流就是运输问题的最优解。

首先,要确认一个初始方案,然后按照最小费用流的负回路算法或最短路径算法进行改进,直到得到最优解为止。

10.2、确认初始方案

作赋权的完全二部图 G=(A,B,E,w)G=(A, B, E, w),其中:

  • A={A1,A2,…,Am}A=\{A_1, A_2, \dots, A_m\},B={B1,B2,…,Bn}B=\{B_1, B_2, \dots, B_n\}
  • E={(Ai,Bj)∣1⩽i⩽m,1⩽j⩽n}E=\{(A_i, B_j)|1\leqslant i \leqslant m, 1\leqslant j \leqslant n\}
  • w:E→Rw: E\rightarrow R

alt text

如果:

{∑j=1nxi,j=ai,1⩽i⩽m∑i=1mxi,j=bj,1⩽j⩽nxi,j⩾0,1⩽i⩽m,1⩽j⩽n\begin{cases} & \sum_{j=1}^n x_{i, j} = a_i, & 1\leqslant i \leqslant m\\ & \sum_{i=1}^m x_{i, j} = b_j, & 1\leqslant j \leqslant n\\ & x_{i, j} \geqslant 0, & 1\leqslant i \leqslant m, 1\leqslant j \leqslant n \end{cases}

则 x={xi,j}x=\{x_{i, j}\} 称为 GG 的一个调运方案。调运方案的总费用为:

w(x)=∑i=1m∑j=1nwi,jxi,jw(x)=\sum_{i=1}^m \sum_{j=1}^n w_{i, j} x_{i, j}

运输问题也可以表述为求使 w(x)w(x) 最小的调运方案。


按照以下方法确定初始调运方案:

  1. 随意选取一条边 (Ai,Bj)(A_i, B_j),然后令 xi,j=min⁡{ai,bj}x_{i, j} = \min\{a_i, b_j\}
    • 如果 ai<bja_i < b_j,则删去 AiA_i 以及所有相关联的边,并令 bj←bj−aib_j \leftarrow b_j - a_i
    • 如果 ai>bja_i > b_j,则删去 BjB_j 以及所有相关联的边,并令 ai←ai−bja_i \leftarrow a_i - b_j
    • 如果 ai=bja_i = b_j,则再选取一条与 AiA_i 或 BjB_j 相关联的边,例如 (Ai,Bk),k≠j(A_i, B_k), k\neq j,然后令 xi,k=0x_{i, k}=0,删去 AiA_i 和 BjB_j 以及所有相关联的边
  2. 重复步骤 1,直到选取了 m+n−2m+n-2 条边为止(注意第一步中保证删去的节点数和选出的边数始终相等),此时还剩下两个顶点 ApA_p 和 BqB_q 和唯一的一条边 (Ap,Bq)(A_p, B_q)且 ap=bqa_p=b_q(产销平衡,第一步中的操作始终保持了这一点)。选取 (Ap,Bq)(A_p, B_q),令 xp,q=apx_{p, q} = a_p。这选中的共 m+n−1m+n-1 条边构成的集合称之为 选中边集,再令除选中边集之外的所有边的 xx 值均为 0。容易验证,这样的 x={xi,j}x=\{x_{i, j}\} 就是一个调运方案。

有命题:按以上方式确定的初始调运方案,选中边集 TT 构成 GG 的一颗生成树。

证明:TT 中有 m+n−1m+n-1 条边,涉及 m+nm+n 个顶点,只需证明 TT 导出的子图是连通的即可。将 TT 的导出子图仍记为 TT。

m+n=2m+n=2 时,m=n=1m=n=1,显然成立。假设对于 m+n=k(k⩽2)m+n=k(k\leqslant 2) 时,结论成立。

现考虑 m+n=k+1m+n=k+1 的情况。按上述方法选取第一条边 (Ai,Bj)(A_i, B_j) 之后,若 ai<bja_i<b_j,GG 删去 AiA_i 以及所有相关联的边,得到子图 G′G',TT 删去 AiA_i 以及 (Ai,Bj)(A_i, B_j) 记作 T′T'。此时 G′G' 仍然是一个有 kk 个顶点的完全二部图,而 T′T' 是 G′G' 中选出来的 k−1k-1 条边构成的子图,由归纳假设,T′T' 是联通的,而 (Ai,Bj)(A_i, B_j) 连接起了 AiA_i 和 T′T',所以 TT 是联通的。若 ai>bja_i>b_j,则类似。

若 ai=bja_i=b_j,不妨设还选中了 (Ai,Bk),k≠j(A_i, B_k), k\neq j,则删去 AiA_i 和 BjB_j 以及所有相关联的边,得到子图 G′G',TT 删去 Ai,BjA_i, B_j 以及 (Ai,Bj),(Ai,Bk)(A_i, B_j), (A_i, B_k) 记作 T′T'。此时 G′G' 仍然是一个有 k−1k-1 个顶点的完全二部图,而 T′T' 是 G′G' 中选出来的 k−2k-2 条边构成的子图,由归纳假设,T′T' 是联通的,而 (Ai,Bj),(Ai,Bk)(A_i, B_j), (A_i, B_k) 连接起了 Ai,BjA_i, B_j 和 T′T',所以 TT 是联通的。证毕。

10.3、改进调运方案

设 T={Air,Bjr∣1⩽r⩽m+n−1}T=\{A_{i_r}, B_{j_r}|1\leqslant r \leqslant m+n-1\}。给每个顶点一个变量,uiu_i 对应 AiA_i,vjv_j 对应 BjB_j。

对每一条选中的边,(Air,Bjr)(A_{i_r}, B_{j_r}),有:

uir+vjr=wir,jr,1⩽r⩽m+n−1u_{i_r} + v_{j_r} = w_{i_r, j_r}, 1\leqslant r \leqslant m+n-1

联立起来,得到一个 m+nm+n 个变量,m+n−1m+n-1 个方程的线性方程组,有无穷多组解。容易找到一组解,这组解中每个顶点的变量值就称为顶点的位势。

记 λi,j=wi,j−(ui+vj)\lambda_{i, j} = w_{i, j} - (u_i + v_j) 为边 (Ai,Bj)(A_i, B_j) 的检验数。显然选中边的检验数均为 0。


调运方案 xx 对应容量-费用网络 NN 上的一个容量为 v0v_0 的可行流,也记作 xx。由于离开 ss 和进入 tt 的边都是饱和边,所以这些边不会出现在 NxN_x 中。因此,NxN_x 中的任意一条回路 CC 都不会经过 ss 和 tt,故 CC 必然是 GG 的一个回路。由于非选中边都是零流边,只能是前向边;选中边既是非饱和边又是非零流边,既可以是前向边又可以是后向边。

由于 TT 是 GG 的一颗生成树,所以 GG 中任意一条非选中边 (Ai,Bj)(A_i, B_j) (非树边)都与 TT 中唯一的一条路径 PP 共同构成一个回路 Ci,jC_{i, j}。记 i0=i,j0=ji_0=i, j_0=j,则 Ci,jC_{i, j} 可以表示为 Ci,j=Ai0Bj0P=Ai0Bj0Ai1Bj1…BjpAi0C_{i, j}=A_{i_0}B_{j_0}P=A_{i_0}B_{j_0}A_{i_1}B_{j_1}\dots B_{j_p}A_{i_0}。

因为从 BB 到 AA 的边都是后向边,权为负。而且除 (Ai,Bj)(A_i, B_j) 的边都是选中边,满足之前的方程,所以 Ci,jC_{i, j} 的权:

w(Ci,j)=wi0,j0−wi1,j0+wi1,j1−wi2,j1+⋯+wip,jp−wi0,jp=wi0,j0−(ui1+vj0)+(ui1+vj1)−(ui2+vj1)+⋯+(uip+vjp)−(ui0+vjp)=wi0,j0−(ui0+vj0)=λi,j\begin{split} w(C_{i, j})&=w_{i_0, j_0}-w_{i_1, j_0}+w_{i_1, j_1}-w_{i_2, j_1}+\dots + w_{i_p, j_p}-w_{i_0, j_p}\\ &=w_{i_0, j_0}-(u_{i_1}+v_{j_0})+(u_{i_1}+v_{j_1})-(u_{i_2}+v_{j_1})+\dots +(u_{i_p}+v_{j_p})-(u_{i_0}+v_{j_p})\\ &=w_{i_0, j_0}-(u_{i_0}+v_{j_0})=\lambda_{i, j} \end{split}

也即 (Ai,Bj)(A_i, B_j) 的检验数,就等于通过它所构成的回路的权,也即通过这个回路调整流量时的单位成本。检验数为负时,通过适当的调整流量,可以降低总费用。


视 (Ai,Bj)(A_i, B_j) 为第一条边,规定 i0=ip+1−i,j0=ji_0=i_{p+1}-i, j_0=j

记调整量:

δ=min⁡{xir+1jr∣0⩽r⩽p}\delta = \min\{x_{i_{r+1}}j_r | 0\leqslant r \leqslant p\}

也即回路中第偶数条边(都是反向边)的最小值。然后做如下调整:

xkl′={xkl+δ,k=ir,l=jrxkl−δ,k=ir+1,l=jrxkl,otherwisex_{kl}'=\begin{cases} x_{kl}+\delta, & k=i_r, l=j_r\\ x_{kl}-\delta, & k=i_{r+1}, l=j_r\\ x_{kl}, & otherwise \end{cases}

也即第奇数条边流量加 δ\delta,第偶数条边流量减 δ\delta,其他边流量不变。

设 ee 是求 deltadelta 时对应的那条 xx 最小的边(若有多条,任取一条),则 x′x' 在 ee 上的值为 0。用 (Ai,Bj)(A_i, B_j) 替换掉 ee,得到新的选中边集 T′T'。

以上操作称之为选中边替换,(Ai,Bj)(A_i, B_j) 称为换入边,ee 称为唤出边。选中边替换,实际上就是找到检验数为负的非选中边,然后通过它所构成的回路调整流量,从而降低总费用。为了保持性质,用 (Ai,Bj)(A_i, B_j) 替换掉流量变为 0 的边。


引理 13:对调运方案 xx 和对应的选中边集 TT 作以 (Ai,Bj)(A_i, B_j) 为换入边的选中边替换得到 x′x' 和 T′T',则 x′x' 仍然是一个调运方案,且 w(x′)=w(x)+λi,j⋅δw(x')=w(x)+\lambda_{i, j}\cdot \delta。

定理 7:一个调运方案时最优的当且仅当所有检验数均不小于 0。

证明:必要性由引理 13 易得。 现证充分性。

设所有的检验数均不小于 0。由定理 5,要证明调运方案最优,只需要证明 GG 中没有权 ww 的负回路。

由前述,所有的基本回路(加入一条非树边构成的回路)均不是负回路。若某个回路 CC 是副回路,则其至少有两条非选中边。任取其中一条非选中边 (Ai,Bj)(A_i, B_j),这条非选中边对应的基本回路记为 Ci,j=X1…AiBj…X2…Y…X1C_{i, j}=X_1 \dots A_iB_j\dots X_2\dots Y \dots X_1,其中 X1,X2X_1, X_2 之间是 CC 和 Ci,jC_{i, j} 的公共部分。C=X1…AiBj…X2…Z…X1C = X_1 \dots A_iB_j \dots X_2 \dots Z \dots X_1。

alt text

令 C′=X2…Z…X1…Y…X2C' = X_2 \dots Z \dots X_1\dots Y \dots X_2,则有 w(Ci,j)+w(C′)=w(C)w(C_{i, j})+w(C')=w(C)。由于 Ci,jC_{i, j} 不是负回路,w(Ci,j)⩾0w(C_{i, j})\geqslant 0,所以 w(C′)⩽w(C)w(C')\leqslant w(C)。因此 C′C' 也是一个负回路,且比 CC 少一条非选中边。重复上述过程,直到得到一个基本回路为止,这与所有基本回路均不是负回路矛盾,证毕。

11、二部图匹配

11.1、定义

简单二部图是一个无向图G=(A,B,E)G=(A, B, E),其中A,BA, B都是顶点集,边集满足∀(i,j)∈E,i∈A,j∈B\forall (i, j)\in E, i\in A, j\in B。

如果存在一个边集M⊂EM\subset E,使得MM中的边两两不相邻(不相邻即指没有公共顶点),则称MM是GG的一个匹配。边数最多的匹配称为最大匹配。当∣A∣=∣B∣=∣M∣|A|=|B|=|M|时,称MM是一个完美匹配。

二部图常用于对两种不同类型的对象进行配对的问题,例如任务分配问题(一个人可以完成某些任务,但是最多只能去完成一个任务,要求分配任务使得尽可能多的人得到任务)。


设MM是二部图GG的匹配,称MM中的边为匹配边,E−ME-M中的边为非匹配边。与匹配边关联的顶点称为饱和点,不与任何一条匹配边关联的顶点称为未饱和点。

GG中,由匹配边和非匹配边交替构成的路径称为交错路径。如果交错路径的起点和终点都是未饱和点,则称这条交错路径是增广交错路径。


引理:设MM是二部图GG的一个匹配,PP是一个增广路径,则:

M′=M⊕EP=M∪EP−M∩EPM'=M\oplus E_P=M\cup E_P-M\cap E_P

是一个匹配,且∣M′∣=∣M∣+1|M'|=|M|+1。其中EPE_P代表路径PP上的边。

例如对以下二部图和其一个匹配(线代表二部图,红线代表匹配):

alt text

则 P=A1B1A2B2A3B4P=A_1B_1A_2B_2A_3B_4 是一条增广交错路径,求异或之后,得到新的匹配(如红线所示):

alt text


由引理易得定理:二部图的匹配是最大匹配当且仅当不存在关于它的增广交错路径。

11.2、匈牙利算法

由上述定理,可以得到匈牙利算法的基本思想:从初始匹配MM开始,不断寻找增广交错路径PP,然后令M←M⊕EPM\leftarrow M\oplus E_P,直到不存在增广交错路径为止。

若(Ai,Bj)∈M(A_i, B_j)\in M,则记match(Ai)=Bj,match(Bj)=Aimatch(A_i)=B_j, match(B_j)=A_i。若AiA_i是未饱和点,则match(Ai)=0match(A_i)=0。


匈牙利算法的完整叙述如下:

  1. 初始化匹配M=∅M=\emptyset
  2. 记AA中已标号未检查的点集为XX,初始化X=∅X=\emptyset
  3. 对于AA中所有未饱和点AiA_i,给其打上标记l(Ai)=0l(A_i)=0,并将其加入XX
  4. 记BB中未标号顶点集为YY,初始化Y=BY=B
  5. 当XX非空时,从XX中取出一个顶点AiA_i,对于YY中AiA_i的每一个邻接点BjB_j(也即BB中与AiA_i邻接且未标记的点),取出,先给BjB_j打上标记l(Bj)=Ail(B_j)=A_i:
    • 如果BjB_j是未饱和点,说明找到了一条增广交错路径,转到步骤 7
    • 如果BjB_j是饱和点,找到BjB_j对应的匹配点AkA_k,将AkA_k加入XX,同时给BjB_j打上标记l(Ak)=Bjl(A_k)=B_j
  6. XX为空,没有找到增广交错路径,返回当前MM,结束算法
  7. 令match(Ai)=Bj,match(Bj)=Aimatch(A_i)=B_j, match(B_j)=A_i
  8. 如果l(Ai)=0l(A_i)=0,说明到达起点,清空所有标记,回到步骤 2;否则,令Bj←l(Ai),Ai←l(Bj)B_j\leftarrow l(A_i), A_i\leftarrow l(B_j),回到步骤 7

定理:匈牙利算法终止时得到的匹配是GG的最大匹配,且算法在O(min⁡{∣A∣,∣E∣}⋅∣E∣)O(\min\{|A|, |E|\}\cdot |E|)时间内终止。

11.3、最大匹配与最大流

二部图的最大匹配问题可以转化为最大流问题。

设G=(A,B,E)G=(A, B, E)是二部图,构造容量网络 NG=(V,E′,c,s,t)N_G = (V, E', c, s, t),其中:

  • V=A∪B∪{s,t}V = A \cup B \cup \{s, t\}
  • E′={(s,Ai)∣Ai∈A}∪{(Ai,Bj)∣(Ai,Bj)∈E}∪{(Bj,t)∣Bj∈B}E' = \{(s, A_i) | A_i \in A\} \cup \{(A_i, B_j) | (A_i, B_j) \in E\} \cup \{(B_j, t) | B_j \in B\}
  • ∀(i,j)∈E′,c(i,j)=1\forall (i, j) \in E', c(i, j) = 1

则 GG 的匹配与 NGN_G 上的 0-1 可行流(边上的流量均为 0 或 1)一一对应,且 GG 的最大匹配的大小等于 NGN_G 上的最大流的流量。因此求最大匹配就等于求最大流。

容易发现,NGN_G 是一个简单容量网络,用 Dinic 算法可以在 O(n1/2m)O(n^{1/2}m) 时间内求解出最大流,从而得到 GG 的最大匹配。