科普文章 - 单位蒙日矩阵乘法

发表于 2023-09-29 分类于科普文章

讲解了单位蒙日矩阵乘法的分治算法的原理和实现

抄写自 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$ 的排列就可以完全确定，我们在这里提前声明好排列和矩阵的对应方式：

对于排列 $p_1,\ldots,p_n$，考虑让（从上到下）第 $i$ 行的黑点处于（从左到右）第 $p_i$ 列，由此生成的单位蒙日矩阵与排列 $p$ 对应。

例如，例图中三个单位蒙日矩阵对应的排列分别是 $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)$ 的位置进行一个分类讨论：

$\delta_{i^+,j^+}\leq0$，那么由于 $\delta$ 的单调性，以 $(i,j)$ 为顶点的四个格子的 $\delta$ 值都是非正的，所以 $C$ 在它们处的取值都是 $C^L$，将无关项约掉之后就得到 $C^M_{i,j}=D^L_{i,j}$；
$\delta_{i^-,j^-}\ge 0$，那么由于 $\delta$ 的单调性，以 $(i,j)$ 为顶点的四个格子的 $\delta$ 值都是非负的，所以 $C$ 在它们处的取值都是 $C^R$，将无关项约掉之后就得到 $C^M_{i,j}=D^R_{i,j}$；
$\delta_{i^+,j^+}>0, \delta_{i^-,j^-}<0$，容易证明相邻格子的 $\delta$ 最多差 $1$，所以这种情况必定是 $\delta_{i^+,j^+}=1, \delta_{i^-,j^-}=-1, \delta_{i^+,j^-}=\delta_{i^-,j^+}=0$，所以我们知道了 $C^L,C^R$ 在这几个位置的差，带入 $C^M$ 即可得到 $C^M_{i,j}=D^{L}_{i,j}+1$，由于 $0\leq C^M_{i,j},D^L_{i,j}\leq 1$，所以只能 $C^M_{i,j}=1$。

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

如果一个蓝点不在红色边线的左上角，那么就保留，否则删除；
如果一个红点不在蓝色边线的右下角，那么就保留，否则删除；
添加所有红蓝边线的交点（记为绿色点）。

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

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

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

代码（loj #184）：

代码

#include <bits/stdc++.h>
using namespace std;
const int MAXN=500010;
vector <int> a,b,c;
int n,d1[MAXN],d2[MAXN],e1[MAXN],e2[MAXN];
int v1[MAXN],v2[MAXN],w1[MAXN],w2[MAXN];
void proc (vector <int> &b) {
    for (int i=1;i<=n;i++) {d1[b[i]]=i;}
    for (int i=1;i<=n;i++) {b[i]=d1[i];}
}
vector <int> mul (vector <int> a,vector <int> b) {
    vector <int> c,a1(1),a2(1),b1(1),b2(1),c1,c2;
    int n=a.size()-1,m=n/2;
    c.resize(n+1);
    if (n==1) {c[1]=1;return c;}
    for (int i=1;i<=n;i++) {
        if (a[i]<=m) {a1.push_back(a[i]);}
        else {a2.push_back(a[i]-m);}
        if (b[i]<=m) {b1.push_back(b[i]);}
        else {b2.push_back(b[i]-m);}
    }
    c1=mul(a1,b1),c2=mul(a2,b2);
    for (int i=0;i<=n;i++) {d1[i]=d2[i]=e1[i]=e2[i]=w1[i]=w2[i]=v1[i]=v2[i]=0;}
    int cnt1=0,cnt2=0;
    for (int i=1;i<=n;i++) {
        if (b[i]<=m) {w1[++cnt1]=i;}
        else {w2[++cnt2]=i;}
    }
    cnt1=cnt2=0;
    for (int i=1;i<=n;i++) {
        d1[i]=d2[i]=0;
        if (a[i]<=m) {d1[i]=w1[c1[++cnt1]];}
        else {d2[i]=w2[c2[++cnt2]];}
    }
    for (int i=1;i<=n;i++) {v1[d1[i]]=v2[d2[i]]=i;}
    for (int i=0,j=n,w=0;i<=n;i++) {
        w+=(d1[i]&&d1[i]>j)+(d2[i]&&d2[i]<=j);
        while (j&&w>=0) {
            w-=(v2[j]&&v2[j]<=i)+(v1[j]&&v1[j]>i),j--;
            if (w<0) {break;}
        }
        e1[i]=j-(w>=0);
    }
    for (int i=0,j=n,w=0;i<=n;i++) {
        w+=(d1[i]&&d1[i]>j)+(d2[i]&&d2[i]<=j);
        while (j&&w-(v2[j]&&v2[j]<=i)-(v1[j]&&v1[j]>i)>0) {
            w-=(v2[j]&&v2[j]<=i)+(v1[j]&&v1[j]>i),j--;
        }
        e2[i]=j+(w<=0);
    }
    for (int i=1;i<=n;i++) {
        if (e1[i-1]+1==e2[i]) {c[i]=e2[i];}
        else {c[i]=d1[i]+d2[i];}
    }
    return c;
}
int main () {
    scanf("%d",&n);
    a.resize(n+1),b.resize(n+1);
    for (int i=1;i<=n;i++) {scanf("%d",&a[i]);}
    for (int i=1;i<=n;i++) {scanf("%d",&b[i]);}
    proc(b);
    c=mul(a,b);
    for (int i=1;i<=n;i++) {printf("%d%c",c[i]," \n"[i==n]);}
    return 0;
}