有向图APSP

r_64 posted @ 2016年12月17日 22:20 in 未分类 with tags 黑科技 , 822 阅读

听说读论文的最好方法是把它翻译一遍。。?(然而我不会严格地翻译。。只会大概地复述一遍。。)这次的论文是《All Pairs Shortest Paths in weighted directed graphs - exact and almost exact algorithms》(Zwick'98)。论文的主题是有向带权图的APSP(All Pair Shortest Path,两两最短路)问题。

背景

两个$n\times n$的矩阵相乘,其naive算法的复杂度是$O(n^3)$。近年来多个更快的算法被提出,最好的bound已经达到了$O(n^{2.376})$。我们用$O(n^\omega)$来表示将两个矩阵相乘的复杂度,那么$\omega<2.376$。很多人猜想$\omega=2+o(1)$,即对任意$\varepsilon>0$,矩阵乘法可以做到$O(n^{2+\varepsilon})$。。当然这个还没被证明。

快速的矩阵乘法是很多最短路算法的基础。比如说无向无权图的最短路可以做到$\tilde{O}(n^{\omega})$。对了,提一下$\tilde{O}$符号。。它表示省略了若干个$\log n$的因子。即,$\tilde{O}(f(n))$可以看做$O(f(n)\text{polylog }n)$。

最短路这个问题其实。。挺难的。据说现在最好的有向实权图(边权是实数)最短路算法也不能做到$O(n^{3-\epsilon})(\forall \epsilon>0)$。论文最后还提了一些open problem,其中就包括“权值绝对值不超过$n$的有向整数权图APSP能否在sub-cubic time内解决?”

再介绍一个符号:$\omega(p,q,r)$表示$n^p\times n^q$的矩阵与$n^q\times n^r$的矩阵相乘的最优复杂度。例如$\omega(1,1,1)=\omega$。$r\le 1$时有一个比较显然的bound:$\omega(1,r,1)\le 2+r(\omega-2)$,就是把问题变成$n^{2(1-r)}$个$n^r\times n^r$的矩阵相乘就好。$\omega(1,r,1)$是有更好的bound的,后面分析复杂度时会讲。

然后介绍Distance Product:两个矩阵$A,B$,不妨令$m=n^r$,$A$的大小是$n\times m$,$B$的大小是$m\times n$。它们的Distance Product被定义为:

$$A\star B=C,c_{ij}=\min_{k=1}^m\{a_{ik}+b_{kj}\}$$

即,考虑一张图分三层,左、右边$n$个点,中间$m$个点,$A$是左边点向中间点连边的权值;$B$是中间向右边连边的权值;那么$A\star B$是左边和右边的最短路权值。

快速计算Distance Product对APSP的计算很有帮助。实数权值的Distance Product应该还没有比$O(n^{2+r})$好到哪里去的算法(注:),不过整数权值的话就可以被规约成矩阵乘法了。。很tricky的一个东西。顺带提一下,FFT大家都知道吧。。计算两个$O(M)$位整数的乘积,复杂度$\tilde{O}(M)$。你们可能已经猜到trick是什么了

这篇论文的results

  • $\tilde{O}(\min\{Mn^{\omega(1,r,1)},n^{2+r})$复杂度的算法,计算一个$n\times n^r$的矩阵和$n^r\times n$的矩阵的Distance Product,要求矩阵中的元素都是绝对值不超过$M$的整数,或者是$+\infty$。
  • 一个$O(n^{2.575})$的随机算法,计算边权范围为$\{-1,0,1\}$的有向带权图的APSP。这么说吧,令$\mu$为$\omega(1,\mu,1)=1+2\mu$的解,则算法的时间复杂度为$\tilde{O}(n^{2+\mu})$。$\mu<0.575$。好像对于小边权的图(而不是简简单单$\{-1,0,1\}$),这个算法也work。对于$M<n^{3-\omega}$来说,算法复杂度是sub-cubic的。
  • 把上面的那个算法derandomization了。也就是说,这篇论文给出了一个$\tilde{O}(n^{2+\mu})$的确定性算法求有向图的APSP。但可惜,这个算法只适用于无权图。听说修改之后可以适用于带权图。
  • 一个复杂度非常优秀的近似算法。它能输出APSP的$1+\epsilon$近似,时间复杂度为$\tilde{O}((n^{\omega}/\epsilon)\log(W/\epsilon))$。$1+\epsilon$近似就是说,如果答案是$a$,输出是$o$,则$a\le o\le a(1+\epsilon)$。图的权值可以是任意实数(近似算法这个条件比较耍赖>_>不过这个复杂度也确实非常非常优秀了)

Distance Product

$\tilde{O}(\min\{Mn^{\omega(1,r,1)},n^{2+r})$计算两个矩阵的Distance Product的算法。那个$n^{2+r}$显然只是来卖萌的,接下来就是用$O(Mn^{\omega(1,r,1)})$计算Distance Product的算法:

algorithm $\mathbf{DIST-PROD}(A,B,M)$
  $a'_{ij}\leftarrow\begin{cases}(m+1)^{M-a_{ij}}&\text{if }|a_{ij}|\le M\\0&\text{otherwise}\end{cases}$
  $b'_{ij}\leftarrow\begin{cases}(m+1)^{M-b_{ij}}&\text{if }|b_{ij}|\le M\\0&\text{otherwise}\end{cases}$
  $C'\leftarrow \mathbf{FAST-PROD}(A',B')$
  $c_{ij}\leftarrow\begin{cases}2M-\lfloor\log_{(m+1)}c'_{ij}\rfloor&\text{if }c'_{ij}>0\\+\infty&\text{otherwise}\end{cases}$

注:上面的算法中$m=n^r$。这个小trick就是将运算挪到指数上,再相加的话主要部分占有绝对的contribution。(例如,考虑$\lfloor\log_5(5^7+5^6+5^6+5^3)\rfloor$,显然会等于$7$,也就是说矩阵乘法中的加法起到了$\max$的作用,乘法起到了$+$的作用)。

算法的正确性就比较显然了。$c'_{ij}=\sum_{k=1}^m(m+1)^{2M-(a_{ik}+b_{kj})}$,其中$a_{ik}=+\infty$或$b_{kj}=+\infty$的$k$不算。那么$2M-\lfloor\log_{m+1}c'_{ij}\rfloor=\min_{k=1}^m\{a_{ik}+b_{kj}\}$。

复杂度的话。。$\mathbf{FAST-PROD}$会调用$O(n^{\omega(1,r,1)})$次整数加法/乘法,而这些运算都是可以$\tilde{O}(M)$做的。故复杂度为$\tilde{O}(Mn^{\omega(1,r,1)})$。

听说将算法改造之后可以返回Witness矩阵。即返回一个$W$使得$c_{ij}=a_{iw_{ij}}+b_{w_{ij}j}$。论文说详情见“full version of this paper”,然而我并没有full version,所以是不是又要留一个小坑了呢?

随机化APSP

大力乱搞不等式证bound要来了。。前方高能预警

先把算法放上来

algorithm $\mathbf{RAND-SHORT-PATH}(D)$
  $F\leftarrow D$;$W\leftarrow 0$
  $M\leftarrow \max\{|d_{ij}|:d_{ij}\ne +\infty\}$
  for $\ell\leftarrow 1$ to $\lceil\log_{3/2}n\rceil$ do
    $s\leftarrow (3/2)^{\ell}$
    $B\leftarrow \mathbf{RAND}(\{1,2,\dots,n\},(9\ln n)/s)$
    $(F',W')\to \mathbf{DIST-PROD}(F[*,B],F[B,*],sM)$
    if $f'_{ij}<f_{ij}$ then $f_{ij}\leftarrow f'_{ij}$;$w_{ij}\leftarrow b_{w'_{ij}}$
  return $(F,W)$

首先声明:算法中的$\mathbf{RAND}(S,p)$返回一个集合,每个元素有$\min(1,p)$的概率在该集合中。$F[*,B]$表示抽出$F$的第$i(i\in B)$列组成的$n\times |B|$的矩阵,$F[B,*]$同理。算法中那个$W$矩阵,储存的就是witness($F_{ij}=F_{iw_{ij}}+F_{w_{ij}j}$)。

这算法在干什么呢?我们考虑一下,如果我们有一个邻接矩阵$D$,怎么通过乱搞来求APSP。我们可以这样乱搞:每一轮选一个点集$B$,把$D_{ij}$与$\min_{k\in B}(D_{ik}+D_{kj})$取较小值——也就是用$S$内的点更新当前求得的APSP。多尝试几轮(讲道理运气非常好的时候,乱选$B$,$\log n$轮就可以得到答案),也许就能得到答案了。

令人惊讶的是。。就算我们的$B$是随机选取的,通过精心选取参数(就是算法中你看到的$(9\ln n)/s$之类的数)我们仍然可以得到一个正确率有严格证明的算法!

令$\delta(i,j)$表示$i$到$j$之间的最短路长度。在算法执行的任何时刻,对任何$i,j\in V$,下面三条总是成立。

  1. $f_{ij}\ge\delta(i,j)$;
  2. $f_{ij}\ne +\infty$,则$f_{ij}=d_{ij}$且$w_{ij}=0$,或者$w_{ij}>0$且$f_{ij}\ge f_{i,w_{ij}}+f_{w_{ij},j}$。
  3. 如果$\delta(i,j)=\delta(i,k)+\delta(k,j)$,且某次迭代的开始有$f_{ik}=\delta(i,k)$,$f_{kj}=\delta(k,j)$,$|f_{ik}|,|f_{kj}|\le sM$,且$k\in B$,则迭代结束之后$f_{ij}=\delta(i,j)$。

1.和2.很容易用归纳法证明,3.的证明如下:迭代结束后$f'_{ij}\le f_{ik}+f_{kj}=\delta(i,j)$,于是$f_{ij}\le \delta(i,j)$,再考虑1.即可证明3.。

更重要,也是更有趣的是下列引理:

令$0\le\ell\le\lfloor\log_{3/2}n\rfloor$,$s_{\ell}=\begin{cases}1&\ell=0\\\lceil 3s_{\ell-1}/2\rceil&\ell>0\end{cases}$,如果$i$和$j$有一条最短路只经过了$s_{\ell}$条边,则第$\ell$次迭代之后,$f_{ij}=\delta(i,j)$的概率非常高。

证明:考虑对$\ell$进行归纳。$\ell=0$的时候,结论显然。我们假设$\ell-1$的时候成立,我们来证明$\ell$的时候也成立。令$p$是$i$到$j$的一条最短路,且$p$有至多$s_{\ell}$条边。如果$p$的边数不超过$s_{\ell-1}$,那么已经成立了,否则我们在$p$上寻找两个点$I,J$(顺序是$i\to I\to J\to j$),使得$i\to I,J\to j$的长度至多是$\lceil s_{\ell}/3\rceil$条边,$I\to J$的长度正好是$\lfloor s_{\ell}/3\rfloor$条边。令$A$为$I\to J$这段路上的所有点组成的集合(包含$I,J$),则$\delta(i,j)=\delta(i,k)+\delta(k,j)$,且$|A|=\lfloor s_{\ell}/2\rfloor+1$,$i\to k,k\to j$的长度都至多为$\lfloor s_{\ell}/3\rfloor+\lceil s_{\ell}/3\rceil\le s_{\ell-1}$(通过讨论可以证明$\le$号的正确性)。并且,$|f_{ik}|,|f_{kj}|\le s_{\ell}M$。根据上面第3条,只要存在$k\in B$的概率比较大,那么结论就是成立的。而$A\cap B=\emptyset$的概率为:

$$P(A\cap B=\emptyset)\le \bigl(1-\frac{9\ln n}{s})^{s/3}\le e^{-3\ln n}=n^{-3}$$

这样,考虑所有点对,我们出错的概率顶多是$n^{-3}\times n^2=n^{-1}$。这就完成了我们的证明。

论文中没有讲。。我算出来正确概率为$(1-\frac{1}{n})^{\lceil\log_{3/2}n\rceil}$。

如果图中没有负权环,我们可以以很高的概率跑出APSP的正确答案,并且通过$W$矩阵构造出最短路径。

下一步是分析复杂度:第$\ell$次迭代的复杂度为$n\times m$的矩阵与$m\times n$的矩阵做Distance Product的复杂度,其中$m=O((n\log n)/s)$,矩阵元素绝对值不超过$sM$。我们令$s=n^{1-r},M=n^t$,那么复杂度为$\tilde{O}(\min\{n^{t+\omega(1,r,1)+1-r},n^{2+r}\})$。注意到$\omega(1,r,1)+1-r$是$r$的减函数,$2+r$是$r$的增函数,那么要得到复杂度瓶颈我们只需要解方程$t+\omega(1,r,1)+1-r=n^{2+r}$就好了。令方程的解为$\mu(t)$,则算法复杂度为$\tilde{O}(n^{2+\mu(t)})$。

特别地,对于权值范围为常数(比如说,$\{-1,0,1\}$)的图,我们可以做到$\tilde{O}(n^{2+\mu})$的复杂度,其中$\mu$是$\omega(1,\mu,1)=1+2\mu$的解。为了说明$\mu<0.575$从而算法是$O(n^{2.575})$的,我们给出$\omega(1,r,1)$的bound:令$\alpha=\sup\{r:\omega(1,r,1)=2+o(1)\}$,则$\alpha\ge 0.294$,且$\alpha\le r\le 1$时$\omega(1,r,1)\le 2+\frac{\omega-2}{1-\alpha}(r-\alpha)+o(1)$。那么$\mu\le\frac{1+\alpha-\alpha\omega}{4-2\alpha-\omega}<0.575$。

还有一点就是,若$\omega=2+o(1)$,则该算法复杂度可以看成$\tilde{O}(n^{2.5})$,这也是本文作者猜测APSP算法的最优复杂度为$\tilde{O}(n^{2.5})$或$\tilde{O}(n^{\omega})$的原因之一。

确定性APSP算法

首先,我们令$\eta(i,j)$为$i$到$j$最短路上的边数,如果有多条最短路取边数最少的那条。

我们考虑$\mathbf{RAND-SHORT-PATH}$为什么需要随机:我们需要找一个很神奇的$B$集合,对$\eta(i,j)$在$(s_{\ell-1},s_{\ell}]$内的所有$(i,j)$,都需要$i\to j$的最短路经过$B$内至少一个点。我们发现这个很像Hitting Set,所以可以用上相关的结论。

从上面的分析看出,我们需要引入一个概念叫做Bridging Set,从而第$\ell$次迭代,$B$只需要是strong $(s/3)$-bridging set就好了。Bridging set是这样定义的:点集$B$被称作$s$-Bridging set,如果$\forall i,j\in V$,若$\eta(i,j)\ge s$,则$\exists k\in B$,使得$\delta(i,j)=\delta(i,k)+\delta(k,j)$。点集$B$被称作strong $s$-Bridging set,如果$\forall i,j\in V$,若$\eta(i,j)\ge s$,则$\exists k\in B$,使得$\delta(i,j)=\delta(i,k)+\delta(k,j)$且$\eta(i,j)=\eta(i,k)+\eta(k,j)$。

那么,要去随机化,我们只需要把$\mathbf{RAND}$函数改成一个确定性的,求一个大小不是很大的strong $(s/3)$-bridging set的算法。我们称它为$\mathbf{FIND-BRIDGE}$。

algorithm $\mathbf{FIND-BRIDGE}(W,s)$
  $\mathcal{C}\leftarrow\emptyset$
  for every $1\le i,j\le n$ do
    $U\leftarrow \mathbf{SUB-PATH}(W,i,j,s)$
    if $|U|\ge s+1$ then $\mathcal{C}\to \mathcal{C}\cup\{U\}$
  return $\mathbf{HITTING-SET}(\mathcal{C})$

$\mathbf{FIND-BRIDGE}$过程接受两个参数:$W,s$。$W$是witness矩阵(用来找最短路径的;我们假设它只能用来找$\eta(i,j)\le s$的点对的最短路,且不保证$\eta(i,j)=\eta(i,w_{ij})+\eta(w_{ij},j)$),$s$表示你要找strong $s$-Bridging set。

$\mathbf{SUB-PATH}$过程根据$W$矩阵,找到至多$s+1$个$i\to j$最短路上的点(如果$W$所代表的那条最短路长度$\le s$则返回所有点;如果$W_{ij}=0$,即$W$没有记录$i\to j$最短路,那么返回$\emptyset$)。显然这个过程可以在$O(s)$内实现。

$\mathbf{HITTING-SET}$过程求$\mathcal{C}$的一个Hitting set。我们随便找篇讲hitting set的近似算法的论文然后套板子就好了。原论文中使用的算法复杂度是$\mathcal{C}$中所有元素大小之和;返回的结果大小至多是$(\ln \Delta)+1$倍最优的fractional hitting-set的大小,其中$\Delta$是出现次数最多的元素的出现次数。考虑每个$V$中的元素选$1/s$的weight,那么$\forall t\in\mathcal{C}$,$|t|\ge s$故这是一个fractional hitting set,且大小顶多为$n/s$。又$\Delta\le n^2$,故我们认为返回的结果大小不超过$n(2\ln n+1)/s$。

$\mathbf{FIND-BRIDGE}$显然是$O(n^2/s)$的。我们考虑只有在$s>n^{0.5}$时使用该算法,否则$B=V$,经过分析,它的复杂度仍然是$\tilde{O}(n^{2+\mu})$的。

注意:这个算法有一个缺陷:$\mathbf{FIND-BRIDGE}$选出来的不一定是strong $s$-bridging set,故这个算法只适用于有向无权图。在这篇论文的full version中将改进这个算法使其适用于带权图。。再说一遍我没有full version。

近似APSP

其实这个算法说来也不玄乎:如果实数Distance Product可以在$O(f(n))$计算出,那么APSP显然可以在$\tilde{O}(f(n))$计算出。所以APSP的难点其实在于将Distance Product提速。考虑近似算法,我们其实只要在做Distance Product的时候,把元素用近似的小整数代替,那么仍然可以$\tilde{O}(n^{\omega})$得到两个矩阵的近似Distance Product,整个复杂度$\tilde{O}(n^{\omega})$也就不那么令人惊讶了。

那么问题的关键浮出水面:一个$(1+\epsilon)$近似的Distance Product算法。我们首先考虑一个将大元素矩阵变成小元素矩阵的算法:

algorithm $\mathbf{SCALE}(A,M,R)$
  $a'_{ij}\leftarrow \begin{cases}\lceil Ra_{ij}/M\rceil&\text{if }0\le a_{ij}\le M\\+\infty&\text{otherwise}\end{cases}$
  return $A'$

算法输入一个矩阵$A$,只考虑其$[0,M]$内的元素,将其转化至$[0,R]$内,并将其他元素变成$+\infty$。有了这个算法,我们就可以考虑我们的Approximate Distance Product的算法:

algorithm $\mathbf{APPROX-DIST-PROD}(A,B,M,R)$
  $C\leftarrow +\infty$
  for $r\leftarrow \lfloor\log_2 R\rfloor$ to $\lceil\log_2 M\rceil$ do
    $A'\leftarrow \mathbf{SCALE}(A,2^r,R)$
    $B'\leftarrow \mathbf{SCALE}(B,2^r,R)$
    $C'\leftarrow \mathbf{DIST-PROD}(A',B',R)$
    $C\leftarrow \min(C,(2^r/R)\cdot C')$
  return $C$

该算法输入两个$n\times n$的矩阵$A,B$,输出一个矩阵$C$,使得:若$\tilde{C}=A\star B$,则$\tilde{c}_{ij}\le c_{ij}\le (1+\frac{4}{R})\tilde{c}_{ij}$。也就是说,$C$是$A\star B$的$1+\frac{4}{R}$近似。下面就来证明这一点。

首先,由于$\mathbf{SCALE}$中我们用的是取上整,所以$\tilde{c}_{ij}\le c_{ij}$。接下来我们证明$c_{ij}\le(1+\frac{4}{R})\tilde{c}_{ij}$。设$c_{ij}=a_{ik}+b_{kj}$。不妨设$a_{ik}\le b_{kj}$,令$2^{s-1}<b_{kj}\le 2^s$,其中$s$是不超过$\log_2 M$的非负整数($b_{kj}=0$时可以认为$s=0$)。$s\le \log_2 R$的话,显然第一轮的时候$c_{ij}$就会被更新为准确值。下面考虑$\log_2 R\le s\le \log_2 M$,考虑$r=s$的这一轮,由$\mathbf{SCALE}$的定义我们有

$$\frac{2^ra'_{ik}}{R}\le a_{ik}+\frac{2^r}{R},\frac{2^rb'_{kj}}{R}\le b_{kj}+\frac{2^r}{R}$$

那么

$$c_{ij}\le \frac{2^ra'_{ik}}{R}+\frac{2^rb'_{kj}}{R}\le a_{ik}+b_{kj}+\frac{2^{r+1}}{R}\le (1+\frac{4}{R})\tilde{c}_{ij}$$

最后一个不等号是因为$2^{r-1}<b_{kj}\le \tilde{c}_{ij}$。这样就完成了正确性证明。

$\mathbf{APPROX-DIST-PROD}$的复杂度为$\tilde{O}(Rn^{\omega}\log M)$。

现在就可以给出最终的$1+\epsilon$-approximate APSP的算法了。

algorithm $\mathbf{APPROX-SHORT-PATH}(D,\epsilon)$
  $F\leftarrow D$
  $M\leftarrow \max\{d_{ij}:d_{ij}\ne+\infty\}$
  $R\leftarrow 4\lceil\log_2 n\rceil/\ln(1+\epsilon)$
  $R\leftarrow 2^{\lceil\log_2 R\rceil}$
  for $\ell\leftarrow 1$ to $\lceil\log_2 n\rceil$ do
    $F'\leftarrow \mathbf{APPROX-DIST-PROD}(F,F,Mn,R)$
    $F\leftarrow \min\{F',F\}$
  return $F$

这个算法就很好理解了:为了使最终算出来的是$1+\epsilon$近似,我们需要调一个$R$,然后做$\log n$次Distance Product就好了。

用数学归纳法可以很简单地证明第$\ell$次迭代算出来的值是原来值的$(1+\frac{4}{R})^{\ell}$近似,其中原来值指的是将$\mathbf{APPROX-DIST-PROD}$换成准确的Distance Product之后算出来的值。这样,最终的精确度为

$$(1+\frac{4}{R})^{\lceil\log_2 n\rceil}\le (1+\frac{\ln(1+\epsilon)}{\lceil\log_2 n\rceil})^{\lceil\log_2 n\rceil}\le e^{\ln(1+\epsilon)}=1+\epsilon$$

(读者们应该都记得$\lim_{n\to\infty}(1+\frac{1}{n})^n=e$)

然后算算复杂度就好了:$R=\tilde{O}(1/\epsilon)$(这里:$\lim_{\epsilon\to 0}\frac{\ln(1+\epsilon)}{\epsilon}=1$),那么复杂度就是$\tilde{O}((n^{\omega}/\epsilon)\log M)$。

如果$\mathbf{DIST-PROD}$返回Witness矩阵,那么$\mathbf{APPROX-SHORT-PATH}$也能返回Witness矩阵。

最后

祝您身体健康。


登录 *


loading captcha image...
(输入验证码)
or Ctrl+Enter