正如我们之前介绍过的, 对于 \((\min, +)\) 的矩阵乘法以及卷积而言, 人类的进展非常缓慢, 目前对于 \(\operatorname{poly}(n)\) 级别的值域, 最快的算法是 Williams 的 \(n^3/\exp \Omega(\sqrt{\log n})\) 复杂度的算法.
接下来我们都忽略复杂度里的 \(n^{o(1)}\) 部分.
对于一些具体的组合问题来说, 它们所调用的 \((\min, +)\) 矩阵乘法有特殊性质. 比如
- Matthew99 利用特殊情况的 \((\min, +)\) 矩阵乘法, 将树的编辑距离问题做到了 \(n^{3-\delta}\) 的时间 (Breaking the Cubic Barrier for (Unweighted) Tree Edit Distance: https://arxiv.org/abs/2106.02026).
- 再比如邓老师, fjzzq2002, V. V. Williams 等人将无向无权图的所有点对最短路的 2-加性近似做到了 \(n^{2.29}\), 改进了几十年前的 \(n^{7/3}\) 做法 (New Additive Approximations for Shortest Paths and Cycles: https://drops.dagstuhl.de/opus/volltexte/2022/16391/pdf/LIPIcs-ICALP-2022-50.pdf)
在组合问题里出现的 \((\min, +)\) 矩阵乘法常常具有如下的性质:
- 对于 \(A*B\), 而言, \(B\) 的每行相邻两项的差是有界的, 也即 \(|B_{i,j} - B_{i,j+1}| \leq O(1)\), 而且两个矩阵的值域都是 \(O(n)\).
于是对于差有界的矩阵的 \((\min, +)\) 乘法的改进就成为了一个重要的问题.
类似的一个问题是值域 \(O(n)\), 单调序列的 \((\min, +)\) 卷积问题. 这两个问题目前最优的结果都是 2022 年段然老师 (然然!) 等人的作品. 分别在矩阵乘法和卷积上做到了 \(n^{(3+\omega)/2}\) 和 \(n^{1.5}\) 的复杂度 (忽略 \(n^{o(1)}\)). 这个问题之前的做法调用了加性组合里的 Balog–Szemerédi–Gowers 定理才做到了 \(n^{1.859}\), 可以说是很大的改进了.
特别地, 他们对矩阵乘法的做法是适用性更加普遍的, 将限制放宽到了行单调矩阵, 也即 \(B_{i,j}\leq B_{i,j+1}\). 注意当 \(|B_{i,j} - B_{i,j+1}|\leq \delta\) 的时候, 将第 \(j\) 行加上 \(\delta j\) 就可以变成单调矩阵, 答案把第 \(j\) 行再减去 \(j\delta\) 就行了.
勾勒
让我们先勾勒出解决这个情况的 \((\min, +)\) 矩阵乘法的大致思路. 我们先假设 \(\omega = 2\), 那么应该得到的复杂度是 \(n^{2.5}\).
我们现在的目的是计算 \(A * B\), 首先考虑 scaling, 我们把两个矩阵的元素除以 \(2\), 得到 \(\lfloor A/2\rfloor\) 和 \(\lfloor B/2\rfloor\). 如果先计算出 \(\tilde C = \lfloor A/2\rfloor * \lfloor B/2\rfloor\), 那么真正的结果 \(C=A*B\) 应该是和 \(2\tilde C\) 很接近的, 答案应该在 \(2\tilde C_{i,j}\) 到 \(2\tilde C_{i,j} + 2\) 之间.
现在我们有了 \(C\) 的一个大致的估计, 想要还原出 \(C\), 我们首先考虑多项式矩阵 \((x^{A_{i,j}})\) 和 \((x^{B_{i,j}})\). 他们相乘得到的矩阵的第 \(i, j\) 项就是所有出现的 \(A_{i,k}+B_{k, j}\). 这个多项式的次数是 \(O(n)\) 的, 这意味着我们直接做有 \(n^3\) 的代价, 不可接受.
但我们只关心其中 \(3\) 项的结果, 所以可以考虑先给多项式的上标 \(\bmod p\), 其中 \(p\) 是一个 \(\Theta(n^{1/2})\) 级别的随机素数. 对于另一组 \(A_{i,k}+B_{k,j}\), 他如果碰到了我们想要探测的下标, 这只有大约 \(1/p\) 的概率, 所以每组 \(i, j\), 我们被影响的错误的项有 \(n/p\) 项, 如果能把他们全都找出来, 把这些贡献减回去, 我们就知道要探测的位置的值到底是多少了.
Comment. 上次见到取整 + 修正误差还是周哥筛.
这个过程最好情况也就是画 \(n^3/p\) 的时间把他们全找出来, 所以我们的这个思路差不多最好就是 \(n^{2.5}\) 的复杂度.
而且问题在于最后一步: 怎么只花 \(n^3/p\) 的时间把它们都找出来呢? 而且我们甚至还没用到单调的条件. 实际上真正的做法过程要稍微复杂一些.
稍弱的算法
让我们先小小的修改一下上面的步骤, 让 \(n^{3-\delta}\) 成为可能.
首先我们不除以 \(2\), 而是除以 \(n^{1/3}\), 这样其实就只有一次 scaling 而不是 log 次. 这样首先值域 \(n^{2/3}\) 的 case 首先显然可以 \(n^{8/3}\) 解决.
接下来我们选取的模数 \(p\) 也是 \(n^{1/3}\) 级别, 这样第二步的多项式矩阵的乘法总大小也是 \(n^{8/3}\).
我们重要的是如何把 \(n^3/p \approx n^{8/3}\) 个错误信息抠出来, 也就是 对于 \(\delta \in [-2, 0]\), 我们要把 \(\lfloor A_{ik}/n^{1/3}\rfloor+\lfloor B_{jk}/n^{1/3}\rfloor = \tilde C_{ij} + \delta + d \cdot p\), 其中 \(d\) 是正整数.
我们考虑先固定一组 \(i, k\), 这个时候我们要比较两个数列, 一个是
另一个是
我们注意: \(\lfloor B_{jk}/n^{1/3}\rfloor\) 是递增数列, 这是我们已知条件. \(\tilde C_{ij}\) 是 \((\min, +)\) 矩阵乘法的结果, 它也是每行递增的. 而它们的值域都是 \(O(n^{2/3})\), 这说明 它们都只有 \(\boldsymbol{O(n^{2/3})}\) 个连续段!
所以接下来对于所有的 \(j\), 我们只需要对每个连续段判一下它们是不是错误信息就好了, 如果是就可以直接把这一段抠出来, 复杂度是 \(n^{8/3} + \#Error\) 的, 也就是期望 \(n^{8/3}\).
Comment. 我觉得这一步还是有点意外的啊, 这个连续段的转化看起来自然, 但是想到原问题是从差有界转换过来, 这个转换居然是这么用的, 就有点奇特了...
真正的算法
真正的算法中, 我们需要用一个更加精妙的方法进行 scaling.
首先取一个随机素数 \(p = \Theta(n^{1/2})\), 我们先经过一些处理, 让 \(a_{i,j} + b_{j, k} \bmod p\) 不进位.
这个东西其实不是很难做到, 首先 \(A\) 矩阵没有什么限制, 而且矩阵整体加一个数也没什么影响, 我们可以拆成 \([0, p/3), (p/3, 2p/3), (2p/3, p)\) 这三个值域的段, 分离出三个矩阵 \(A_0, A_1, A_2\), 每个数只塞到一个地方, 剩下的地方取一个大于 \(n\) 且 \(\bmod p\) 正确的数. 只要算出它们各自和 \(B\) 的乘积就能凑出答案了. 对于 \(B\) 矩阵来说我们如法炮制, 但由于有单调性这个限制, 我们不能填 \(+\infty\), 而是基于单调性往上取整到合适的数. 注意所有位置的贡献最后会取 \(\min\), 所以这些往上取整的数是不会影响答案的.
这样一来, 考虑 \(\lfloor A/p\rfloor\) 和 \(\lfloor B/p\rfloor\) 的乘积的时候, 我们得到的 \(\tilde C=\lfloor A/p\rfloor * \lfloor B/p\rfloor\) 就恰好是 \(C=A*B\) 的除以 \(p\) 下取整的结果了.
现在重新审查这个问题, 我们相当于分解成了一个二维问题: \(A,B\) 的值域是一个二元组 \((\lfloor A_{ij}/p\rfloor, A_{ij}\bmod p)\), 由于没有进位, 我们的比较就是直接按照字典序.
接下来的叙述里, 我们把 \(A, B\) 都视作二元组.
我们首先算出关于第一元的最优解, 这只需要 \(n^{2.5}\) 的时间.
对于第二元, 我们接下来考虑 scaling, 先算出把这一维 \(/2\) 的结果, 然后考虑乘 \(2\) 之后的情况.
我们先只把第二元挂在矩阵上, 这是 \(n^{0.5}\) 次的多项式矩阵, 乘法是 \(n^{2.5}\) 的. 然而, 这些乘法的结果包含了第一元没有达到最优的那些 \(A_{i,k}+B_{k,j}\), 我们需要把它们去掉才能提取出答案.
由于我们已经计算出了 \(C\) 的第一维, 我们考虑枚举 \(i, k\), 那么对于 \(j\) 来说, 由于 \(B_{k,j}\) 是单调且值域 \(n^{0.5}\) 的, \(C_{i,j}\) 也是如此, 总共我们能够从矩阵里抠出 \(n^{2.5}\) 个连续段, 这些连续段如果全都去掉它们的贡献, 我们就算完了.
但这些贡献依然是难以处理的, 所以需要考虑一个更加迂回的策略: 我们要考虑当前 \(B_{i,j}\) 这个二元组的连续段, 但它的数量和第二维的大小也有关系了. 但所幸我们要处理的是连续段里会影响到我们查询的多项式次数位置的下标是 \(O(1)\) 范围的, 我们只关心那部分, 所以在 scaling 的过程中我们把之前处理的连续段下传, 然后只保留值域在对应的 \(A_{ik}+B_{kj} = C_{ij} + O(1)\), 且第一维加起来不等 的那些段.
首先看看有了这些段之后我们可以做什么. 这些段的总长可以是很长的, 但庆幸的是, 我们要消除一段 \(A_{i,k}+ B_{k,j}, j\in [l, r]\) 的贡献, 只需要在乘法出来的那个矩阵的 \(i, [l, r]\) 这段打个标记就行了. 所以代价就是连续段的总数.
递归上的时候看起来连续段会越来越多, 但注意到我们始终有个剪枝是要求 \(A_{ik}+B_{kj} = C_{ij} + O(1)\).
从值域的角度论述, 连续段将会有 \(n^3 / 2^l\) 个, 但在现在 \(n/2^l\) 的值域来说, 可以看做随机选取了一个 \(p/2^l\) 的模数, 模完之后的部分相等了, 这样的连续段只有 \((p/2^l)^{-1}\) 的概率保留, 所以连续段的数量还是 \(n^3/p\).
这样一来, 复杂度就是 \(n^{2.5}\) 了.
上面对于期望的叙述是可以严格表述的, 也就是因为一个数只有 \(O(\log n)\) 个质因子.
冷静计算复杂度
刚刚我们一直都假设了 \(\omega = 2\).
对于真正的复杂度, 我们看看每一步 \(\lfloor A/p\rfloor * \lfloor B/p\rfloor\) 到底是怎么做的. 设 \(p=n^{\alpha}\), 对于第一步来说, 用矩阵乘法会得到 \(n^{\omega + 1 - \alpha}\), 但实际上这里还是不需要, 因为我们还是可以根据连续段的数量, 用个数据结构来维护, 做到 \(n^{3-\alpha}\).
矩阵乘法的部分都是 \(n^{\omega+\alpha}\) 的形式, 其他部分用到连续段的依然是 \(n^{3-\alpha}\).
平衡一下就是 \(\alpha=(3-\omega)/2\), 也就是复杂度 \(n^{(3+\omega)/2}\).
单调序列的卷积
容易把上面的做法稍作修改, 得到 \(n^{1.5}\) 的 \((\min, +)\) 卷积. 读者可自行验证.