单位蒙日矩阵乘法

· · 个人记录

代码：https://loj.ac/s/1477463

怎么有人隔三岔五模拟赛里出自己论文题啊？

抄写自 IOI 2022 论文集《浅谈一类蒙日矩阵》。

这篇论文前面是一些基础概念，然后引出了蒙日矩阵的 min+ 乘法。一般蒙日矩阵的 O(n^2) min+ 乘法（SMAWK）之前 Itst 在 APIO 就讲过了，而且和本文也没啥关系，就不说了，我们快进到简单单位蒙日矩阵乘法这里。

前置定义

我们用 i^+,i^- 分别表示 i+0.5,i-0.5，虽然感觉不太习惯，但是这么定义确实使得之后的定义很形象。

定义 1：对于 (n+1)\times (m+1) 矩阵 A，定义一个 n\times m 矩阵 B，满足 B_{i,j}=A_{i^+,j^-}+A_{i^-,j^+}-A_{i^+,j^+}-A_{i^-,j^-}，记为 B=A^M,\ A=B^D，称 B 是 A 的密度矩阵，A 是 B 的分布矩阵。

例如若 A 的行标集合是 0,1,\ldots,n，那么 B 的行标集合实际上是 0.5,1.5,\ldots,n-0.5。

我们可以形象地理解：将 A 中的数记在一个 (n+1)\times (m+1) 表格的每个格子中，B 中的数就记在表格里非边界的格点上，这个数其实是表格从左下到右上的差分，也就是说格子里的数等于其左下角所有格点处数之和。

定义 2（蒙日矩阵）：如果矩阵 A 满足 A^M 的所有元素非负，则称其为蒙日矩阵。如果 A 的第一行和最后一列全零，则称其为简单蒙日矩阵。

定义 3（置换矩阵）：每行每列恰好有一个 1，其他位置全零的矩阵。

定义 4（单位蒙日矩阵）：如果矩阵 A 满足 A^M 是置换矩阵，则称其为单位蒙日矩阵。

定义 5（亚单位蒙日矩阵）：如果矩阵 A 满足 A^M 的每行每列有最多一个 1，其余位置都是 0，则称其为亚单位蒙日矩阵。

定义 6（(\min,+) 乘法）：有 n\times m 矩阵 A 和 m\times p 矩阵 B，如果 n\times p 矩阵 C 满足 \forall i,j,\ C_{i,j}=\min(A_{i,k}+B_{k,j})，则称 C 是 A,B 的 (\min,+) 乘法的乘积，记为 C=A\times B，论文中叫做距离乘法。

首先来感性理解一下单位蒙日矩阵，按照上面说的形象认识，单位蒙日矩阵实际上就是在一个 (n+1)\times (n+1) 网格中，非边界的格点中出现了 n 个黑点，并且每行每列恰好有一个黑点，然后我们在格子中写数，每个格子中的数就是其左下角黑点的个数，下面就展示三个单位蒙日矩阵，黑点的位置也已经标出：

我们还发现，左边两个矩阵进行 (\min,+) 乘法后恰好就等于右边的矩阵，然后我们开始思考，是否任意两个同大小的单位蒙日矩阵乘积还是单位蒙日矩阵呢？答案是肯定的，证明留在附录。

大小为 n 的单位蒙日矩阵只需要一个 1,\ldots,n 的排列就可以完全确定，我们在这里提前声明好排列和矩阵的对应方式：

例如，例图中三个单位蒙日矩阵对应的排列分别是 3\ 1\ 2\ 4,\ 2\ 1\ 4\ 3,\ 4\ 2\ 1\ 3。

下面要给出一种复杂度为 O(n\log n) 的计算两个大小为 n 的单位蒙日矩阵乘积的算法，其输入输出都是与单位蒙日矩阵对应的排列（否则空间就要 O(n^2)）。

算法框架

先说一下，算法中有些部分可能有点绕，如果理解不清楚，建议自己写一下涉及的每个矩阵分别是几乘几的，是置换矩阵，单位蒙日矩阵还是亚单位蒙日矩阵。

计算简单单位蒙日矩阵 A^D,B^D 的乘积（也就是 A,B 分别是置换矩阵），基本想法是分治，我们设答案为 C=A^D\times B^D，那么我们知道：

C_{i,j}=\min_{k\in [1,n]}(A^D_{i,k}+B^D_{k,j})

现在我们将 k 的范围分成 [1,m] 和 (m,n] 两个部分，其中 m 取 \dfrac{n}{2}，然后：

C^{L}_{i,j}=\min_{k\in [1,m]}(A^D_{i,k}+B^D_{k,j}),\ \ C^{R}_{i,j}=\min_{k\in (m,n]}(A^D_{i,k}+B^D_{k,j}) \\ C_{i,j}=\min(C^{L}_{i,j},C^{R}_{i,j})

但是这里其实 C^L,C^R 没有太好的性质，不太好直接算，但是至少我们看到了可以将 A,B 分别分成式子中 k\in [1,m] 和 k\in (m,n] 这两个部分，也就是 A 可以划分成左右两半，左边叫 A^L，右边叫 A^R。而 B 可以划分成上下两半，上边叫 B^L，下边叫 B^R，实际上 C^L 就是 A^{DL} 和 B^{DL} 的乘积，只不过因为 A^{DL} 和 B^{DL} 性质不好，所以才没法直接算。

不过我们换个角度，A^{DL},B^{DL} 性质不好，但是 A^{LD},B^{LD}（分别是 n\times m 和 m\times n 矩阵）性质挺好的，因为它们分别是简单亚单位蒙日矩阵（A^L,B^L 每行每列只有至多一个 1），而删除（A^L 中的）全零行和（B^L 中的）全零列之后（都变成了 m\times m 矩阵）就变成了简单单位蒙日矩阵，可以记一个 D^{LD}=A^{LD}B^{LD}，即 D^L=(A^{LD}B^{LD})^M，同理 D^R，这个形式就和原问题一模一样，直接递归就可以算出来（A^L,B^L,A^R,B^R 里 1 的个数只有原来一半），这里 D^L,D^R 是 m\times m 和 (n-m)\times (n-m) 的置换矩阵。

然后我们再融合 D^L,D^R，例如 D^L 中的 1 就按照原本 A^L 中 1 的行和 B^L 中 1 的列嵌入回 n\times n 的大矩阵中，D^R 同理，这样就拼成了一个新的 n\times n 的简单单位蒙日矩阵的密度矩阵，记为 D。

上面说的有点绕，看一个图例，图中红色的表示抽出的 A^L,B^L 以及递归算出的 D^L，蓝色是 A^R,B^R,D^R，最后嵌入得到 D。

但这个 D 不是我们最终答案，毕竟首先我们把 A^{DL} 换成了 A^{LD}，而且只是对置换矩阵直接合并，怎么想都不是很有道理，所以还需要经过一些神秘变换才能得到最终答案 C（或者说对应的排列 C^M）。

神秘变换

下面来推一下 C^L,C^R 的式子，进而合并出 C：

\begin{aligned} C^L_{i,j} & =\min_{k\in [1,m]}(A^D_{i,k}+B^D_{k,j}) \\ & =\min_{k\in [1,m]}(A^{LD}_{i,k}+B^{LD}_{k,j}+B^{RD}_{m+1,j})\qquad\text{（前缀和拆成上下两部分）}\\ & =\min_{k\in [1,m]}(A^{LD}_{i,k}+B^{LD}_{k,j})+B^{RD}_{m+1,j}\\ & =D^{LD}_{i,j}+D^{RD}_{0,j} \end{aligned}

在最后一步我们指出 B^{RD}_{m+1,j}=D^{RD}_{0,j}，这是因为 D^{RD}_{0,j}=\min_{k\in (m,n]}(A^{RD}_{0,k}+B^{RD}_{k,j})，而 A^{RD}_{0,k}=k-m-1，所以 \min 中的式子关于 k 递增（实际上不是很困难，可以画图理解一下），所以 k=m+1 时取最小值。

同理我们有：

C^{R}_{i,j}=D^{RD}_{i,j}+D^{LD}_{i,n}

然后我们就可以开始考察 C_{i,j}=\min(C^L_{i,j},C^R_{i,j})，考虑 \delta_{i,j}=C^L_{i,j}-C^R_{i,j}，当 \delta_{i,j}\ge 0 时取 C_{i,j}=C^R_{i,j}，否则取 C_{i,j}=C^L_{i,j}，然后推一下式子：

\begin{aligned} \delta_{i,j} & =C^L_{i,j}-C^R_{i,j}\\ & =(D^{LD}_{i,j}+D^{RD}_{0,j})-(D^{RD}_{i,j}+D^{LD}_{i,n})\\ & =(D^{RD}_{0,j}-D^{RD}_{i,j})-(D^{LD}_{i,n}-D^{LD}_{i,j}) \end{aligned}

我们已经说过，D^L,D^R 是一个 n\times n 置换矩阵 D 的两个部分，还是用最开始的黑点模型，考虑 D 中，把来自 D^L,D^R 中的点分别叫做红点和蓝点，那么 D^{LD}_{i,j},D^{RD}_{i,j} 就分别是 (i,j) 左下角的红点个数和蓝点个数，于是 D^{RD}_{0,j}-D^{RD}_{i,j} 就是 (i,j) 左上角的蓝点个数，D^{LD}_{i,n}-D^{LD}_{i,j} 就是 (i,j) 右下角的红点个数，那么 \delta_{i,j} 就是 (i,j) 左上角的蓝点个数减去右下角的红点个数，将它记在每个格子中，如下图 (a)：

由于我们需要决定每个位置取 C^L 还是 C^R，所以实际上我们关心 \delta_{i,j} 的正负，我们发现 \delta_{i,j} 从左上到右下递增，因此我们可以划出两条左下到右上的折线，分出 \delta_{i,j}<0,\ \delta_{i,j}=0,\ \delta_{i,j}>0 的三个部分，如图 (b)，红线左上是 >0，蓝线右下是 <0，中间是 =0。

由于我们定义的是 C=A^D\times B^D，不过最后实际上要输出一个 C^M 对应的排列，所以我们要考虑的就是 C^M 中哪些位置为 1，而 C^M_{i,j}=C_{i^+,j^-}+C_{i^-,j^+}-C_{i^+,j^+}-C_{i^-,j^-}，因此我们对格点 (i,j) 的位置进行一个分类讨论：

然后我们形象地解释一下上面三种分类：

于是对 (b) 图进行上述操作，就得到了 (c) 图，这就是 C^M，也就是我们的答案了（根据每行哪一列有点得到答案排列）。

整个过程由于每次分成两半递归，其余部分都是线性，所以复杂度为 O(n\log n)。

原论文里的流程图（由于集训队论文没电子版，不能抄 dengyaotriangle 的图了！）：

扩展

本质上我们定义了一种排列上的二元运算，输入两个排列输出它们对应的简单单位蒙日矩阵乘积对应的排列，这种运算还有一些内涵和性质。

后面再写。

附录&补充

后面再写。