有向图APSP

听说读论文的最好方法是把它翻译一遍。。?(然而我不会严格地翻译。。只会大概地复述一遍。。)这次的论文是《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的算法：

注：上面的算法中$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要来了。。前方高能预警

先把算法放上来

首先声明：算法中的$\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}$。

$\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算法。我们首先考虑一个将大元素矩阵变成小元素矩阵的算法：

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

该算法输入两个$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的算法了。

这个算法就很好理解了：为了使最终算出来的是$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矩阵。

最后

祝您身体健康。