[论文阅读报告] All pairs shortest paths using bridging sets and rectangular matrix multiplication, FOCS '98

本篇文章介绍整数边权下 \((\min, +)\) 矩阵乘、APSP 等问题的一些做法。若每个元素的权值在 \([-M, M] \cap \mathbb Z\) 中，\(n \times n^r\) 和 \(n^r \times n\) 的 \((\min, +)\) 矩阵乘可做到 \(\tilde O(Mn^{\omega(r)})\)；有向图 APSP 可做到 \(\tilde O(n^{2 + \mu(t)})\)，其中 \(M = n^t, \mu(t)\) 为 \(\omega(r) = 1 + 2r - t\) 的解，当 \(M = O(1)\) 时为 \(\tilde O(n^{2.575})\)；\((1 + \epsilon)\) 近似非负实数 \((\min, +)\) 矩阵乘可以做到 \(\tilde O(n^{\omega(r)}\log W / \epsilon)\)，其中 \(W\) 为最大权值和最小权值的比。

如果我们想要精确计算实数矩阵上的 \((\min, +)\) 矩阵乘法，那么目前不存在 \(O(n^{3 - \delta})\) 的算法，这被称为 APSP hypothesis。但是如果每个元素都是 \([-M, M]\) 中的整数，则我们有很简单的 \(O(n^{3 - \delta})\) 的做法，可以做到 \(\tilde O\left(\min\left(n^{2 + r}, M n^{\omega(r)}\right)\right)\)，其中 \(\omega(r)\) 为 \(n \times n^r\) 和 \(n^r \times n\) 矩阵相乘的快速矩阵乘法指数。

令 \(m = n^r\)。我们让 \(a'_{i, j} = (m + 1)^{M - a_{i, j}}, b'_{i, j} = (m + 1)^{M - b_{i, j}}\)，计算 \(C' = A' \times B'\)，令 \(c_{i, j} = 2M - \operatorname{msb}(c_{i, j}')\)，其中 \(\operatorname{msb}\) 表示 \(c_{i, j}'\) 最高位的位置，也即 \(\lceil \log_{m + 1} c'_{i, j} \rceil\)。不难验证这个做法是正确的，并且实际上计算出了每一种 \(a + b\) 的个数。这样每一次代数运算是 \(\tilde O(M)\) 的，需要的代数运算次数是 \(O(n^{\omega(r)})\) 的，因此复杂度为 \(\tilde O(M n^{\omega(r)})\)。如果 \(M\) 太大，也可以直接暴力计算，复杂度为 \(\tilde O(n^{2 + r})\)。

我们能在 \(\tilde O(M n^{\omega(r)})\) 计算每个位置 \((i, j)\) 中每一种 \(a + b\) 的个数，似乎是不寻常的，因为矩阵乘法加速的原理在于压缩信息，但我们似乎保留了全部的信息。这是由于 \(\tilde O(M n^{\omega(r)})\) 实际上并不是一个多项式复杂度，输入规模是 \(O(n^2 \log M)\) 的。这就类似于桶排序，原本长度为 \(n^r\) 的信息变为了长度为 \(M\) 的，因此保留信息只需要多付出 \(\tilde O(M)\) 倍的时间，而 \(O(n^{\omega(r)})\) 是不变的。

如果我们同时考虑复杂度中的 \(M\) 和 \(n\)，那么 \((\min, +)\) 矩阵乘和 APSP 问题不是等价的。如果使用原先的归约，即使 APSP 问题中每条边边权在 \([-M, M] \cap \mathbb Z\) 中，在 \(\log_2 n\) 轮 \((\min, +)\) 矩阵乘中，权值上限将不断增加，最终达到 \(nM\)，这样做 APSP 还不如实数边权 Dijkstra 快。

不过 APSP 有其图论性质。一条性质是，当最短路用到的边数很多时，它可选择的松弛点也会变多。我们进行 \(\log_{3/2} n\) 轮 \((\min, +)\) 矩阵乘，第 \(i\) 次只考虑松弛那些走过的边的条数不超过 \((3/2)^i\) 的最短路。可以观察到，可用的松弛点的个数是不小于 \(1/3\) 其长度的。当最短路长度很长时，我们可以随机从 \(V\) 中选取一部分点作为可能的松弛点。

假设我们在做第 \(i\) 轮 \((\min, +)\) 矩阵乘，当前的矩阵是 \(F\)，令 \(s = (3/2)^i\)，将 \(V\) 中每个点以 \(\min(1, (9 \ln n) / s)\) 的概率选进集合 \(B\) 中，计算 \(F' = F[*, B] \times F[B, *]\)，将 \(F'\) 中超过 \(sM\) 的元素设为 \(+\infty\)。考虑那些长度处于 \((2s/3, s]\) 的最短路 \((u, v)\)，使用归纳法，如果所有长度 \(\leq 2s/3\) 的最短路已经在 \(F\) 中，此时令 \(P_{uv}\) 中 \(u\) 向后 \(2s/3\) 个点的位置为 \(y\)，\(v\) 向前 \(2s/3\) 个点的位置为 \(x\)，则 \((x, y)\) 间隔了 \(2 \cdot 2s/3 - s = s/3\) 个点，而对所有 \(z \in P_{xy}\)，\((u, z), (z, v)\) 的距离不超过 \(2s/3\)，因此已经在 \(F\) 中，只要有至少一个 \(z \in B\)，\((u, v)\) 的最短路就在 \(F'\) 中。这样的概率是

\[\left(1 - \frac{9 \ln n}s\right)^{s / 3} \leq \exp(-3 \ln n) = n^{-3} \]

使用 union bound，所有长度在 \((2s/3, s]\) 中的 \((u, v)\) 的最短路 \(P_{(2s/3, s]}\) 都在 \(F'\) 中的概率不小于 \(1 - |P_{(2s/3, s]}|n^{-3}\)。将每一轮的概率一起使用 union bound，成功概率不小于 \(1 - P_{[0, nM]} \cdot n^{-3} = 1 - n^{-1}\)。

设 \(s = n^{1 - r}, M = n^t\)，则这一轮的复杂度为 \(\tilde O(\min(sM n^{\omega(r)}, n^{2 + r})) = \tilde O(\min(n^{t + \omega(r) + (1 - r)}, n^{2 + r}))\) 其关于 \(s\) 的最大值在 \(\omega(r) = 1 + 2r - t\) 取到。当 \(M = O(1)\) 时，查表可得，复杂度为 \(\tilde O(n^{2.575})\)（2002，论文发布时）\(\tilde O(n^{2.528})\)（2023，最新结果）。

对于 \(r = 1\) 的情况，我们观察这个式子，发现当 \(M = O(n^{3 - \omega - \delta})\) 时，这个算法是 \(\tilde O(n^{3 - \delta})\) 的。由于 \(\omega \geq 2\)，当 \(M = O(n)\) 时，我们没法从中得到 \(O(n^{3 - \delta})\) 的做法，因此他不违反 APSP 经典假设。

论文继续介绍了构造最短路与去随机化的方法，这些部分较为独立，这里不再赘述。

对于非负实数近似 \((\min, +)\) 矩阵乘，我们可以利用加法的特性给出一个十分经典的近似算法，取 \(R \ll M\)，将权值上取整，即 \(a'_{i, j} = \left\lceil \frac{Ra_{i, j}} M\right\rceil\)，进行同样的计算后再转换回来。

对于不同大小的 \(a_{i, k} + b_{k, j}\)，我们不能把他们都用 \(M\) 作为分母取整，因为在计算近似比时我们使用的是结果的值作为分母而非 \(M\)，这样会导致太小的 \(a + b\) 近似比太差，我们需要枚举 \(r \in [\lfloor \log_2 R \rfloor, \lceil \log_2 M \rceil] \cap \mathbb Z\)，用 \(2^r\) 做分母。当 \(2^{r - 1} < \max(a_{i, j}, b_{j, k}) \leq 2^r\) 时，有

\[a_{i, j} \leq \frac{2^r a'_{i, j}}R < a_{i, j} + \frac{2^r}R,\ b_{j, k} \leq \frac{2^r b'_{j, k}}R < b_{j, k} + \frac{2^r}R \]

\[c_{i, j} \leq c'_{i, j} \leq \frac{2^r a'_{i, j}}R + \frac{2^r b'_{j, k}}R < a_{i, k} + b_{k, j} + \frac{2^{r + 1}}R = c_{i, j} + \frac{2^{r + 1}}R \leq \left(1 + \frac 4R\right) c_{i, j} \]

令 \(R = \frac 4{\epsilon}\) 即可做到 \(1 + \epsilon\) 近似。复杂度为 \(\tilde O(n^{\omega}\log M / \epsilon)\)。对于 APSP，由于我们要做 \(\log n\) 轮，每一轮的误差会相乘，即 \(\left(1 + \frac 4R\right)^{\lceil \log n \rceil}\)，令 \(R = O\left(\frac{\log n}{\ln(1 + \epsilon)}\right) = \tilde O\left(1/\epsilon\right)\) 即可。

我们看到，将 \(M\) 变为 \(4 \log M / \epsilon\) 利用的是不同大小的 \(a + b\) 的采样密度不同，我们总是把 \((2^{r - 1}, 2^r]\) 分成 \(R\) 份，因此采样密度是指数级增长的。在无权图上这一点会更加容易观察到：我们只用计算 \(01\) 矩阵乘法 \(A^{\lfloor (1 + \epsilon)^i \rfloor}, A^{\lceil (1 + \epsilon)^{i + 1} \rceil}\) 看 \((u, v)\) 的值是否从 \(0\) 变为 \(1\) 便可以知道是否有 \(\delta(u, v) \in (\lfloor (1 + \epsilon)^i \rfloor, \lceil (1 + \epsilon)^{i + 1} \rceil]\)。复杂度为 \(O(n^\omega \log_{1 + \epsilon} n) = \tilde O(n^{\omega} / \epsilon)\)。

如果原先权值是非负实数，我们可以将所有的权值放缩，使得最小的权值为 \(1\)（实际上由于边界问题，\(2\) 更好）。这样便囊括在了 \(2^r\) 的枚举中。