题解：P1753 矩阵链排序问题
Chegada_
·
·
题解
题解：P1753 矩阵链排序问题
可能是缘分吧，我在随机跳题的时候刚好跳到了这道题，
看一眼数据范围，n\le 2\times 10^6。可是我只会 O(n^3) 的做法，于是我迅速的打开这道题的题解区，准备借鉴一下大佬的代码，可是空空如也，于是我下定了决心。
ps：因为网上关于 Hu-Shing 算法的资料非常之少，作者只能自己摸索，代码可能不是很严谨，没有实现的很完美，见谅。
## Descreption
- 给定 $n$ 个矩阵按顺序连乘 $M_1 M_2\cdots M_n$，第 $i$ 个矩阵尺寸为 $w_{i-1}\times w_i$（输入给出 $w_0,w_1,\dots,w_n$ 共 $n+1$ 个数）。
- 一个 $a\times b$ 的矩阵乘一个 $b\times c$ 的矩阵，代价为 $a\cdot b\cdot c$。矩阵乘法满足结合律，请通过合理加括号，最小化整条链乘完的总代价。
- 数据范围：$1\le n\le 2\times 10^6$，$1\le w\le 10^4$。其中 $30\%$ 满足 $n\le 500$，另 $30\%$ 满足 $n\le 2\times 10^5$。
因为 $w$ 最大 $10^4$、三角形代价最大 $(10^4)^3=10^{12}$、三角形数量约 $2\times 10^6$，答案最大可到 $2\times 10^{18}$，所以中间运算要用 `__int128` 实现较为稳妥，这一点后面会反复用到。
## Preparations
**先上前置芝士**：
- [区间DP](https://www.luogu.com.cn/problem/P1775)。
- [四边形不等式/Knuth 优化概念](https://oiwiki.moe/dp/opt/quadrangle/)。
- [四边形不等式/Knuth 优化例题](https://www.luogu.com.cn/problem/P4767)。
- [双向链表](https://www.luogu.com.cn/problem/B4324)。
- [可并堆](https://www.luogu.com.cn/problem/P3377)。
如果你准备好了，那么跟着我的思路开始吧！
## 区间 DP（30pts）
这是几乎所有人的第一反应，也是《算法导论》的教科书例题。
记 $f(l,r)$ 为把 $M_l,M_{l+1}\cdots M_r$ 乘起来的最小代价。枚举最后一次乘法发生在哪里（即最外层括号的分界点 $k$），左半 $M_l\cdots M_k$ 得到一个 $w_{l-1}\times w_k$ 的矩阵，右半 $M_{k+1}\cdots M_r$ 得到一个 $w_k\times w_r$ 的矩阵，二者相乘再花 $w_{l-1}\cdot w_k\cdot w_r$：
$$
f(l,r)=\min_{l\le k<r}\Big\{f(l,k)+f(k+1,r)+w_{l-1}\,w_k\,w_r\Big\}, f(i,i)=0
$$
答案就是 $f(1,n)$。按区间长度从小到大转移即可。
```cpp
#include<bits/stdc++.h>
#define int long long
using namespace std;
const int N=510;
int n,i,j,k,l,r,len;
int w[N],f[N][N];
signed main()
{
ios::sync_with_stdio(0);
cin.tie(0);
cin>>n;
for(i=0;i<=n;i++)
cin>>w[i];
for(len=2;len<=n;len++)
{
for(l=1;l+len-1<=n;l++)
{
r=l+len-1;
f[l][r]=4e18;
for(k=l;k<r;k++)
{
f[l][r]=min(f[l][r],f[l][k]+f[k+1][r]+w[l-1]*w[k]*w[r]);
}
}
}
cout<<f[1][n];
}
```
时间复杂度 $O(n^3)$，空间 $O(n^2)$。$n\le 500$ 稳过，但是再往上就不行了，所以这份代码只能吃到那 $30$ 的分，大概**普及**难度吧。
## 一个想当然的坑：四边形不等式在这里用不了
很多人（包括当时的我）看到**区间 DP+求最优分割点**，第一反应是上**四边形不等式/Knuth 优化**，把这个一套，$O(n^3)$ 变 $O(n^2)$。
我兴冲冲写完，样例过了，然后开始心虚 —— 于是写了个暴力对拍，专门检验决策点是否单调，即是否恒有
$$
\mathrm{opt}(l,r-1)\ \le\ \mathrm{opt}(l,r)\ \le\ \mathrm{opt}(l+1,r).
$$
结果**很快就被打脸**。这里给一个具体反例（我从随机数据里抓出来的最小之一），相邻五个维度是
$$
\dots,\ w=2,\ 25,\ 18,\ 23,\ 9,\ \dots
$$
在对应区间上算出来：$\mathrm{opt}(l,r)=6$，而 $\mathrm{opt}(l+1,r)=5$。也就是**区间右端不动、左端右移一格，最优分割点反而往左跳了**，这明显违反单调性，而且不是取值相同，而是货真价实的反例。
原因我找了很久：Knuth 优化要求转移里那项**合并代价 $w(l,r)$ 只与区间 $(l,r)$ 有关**。但矩阵链里的合并代价是 $w_{l-1},\,w_k,\,w_r$，它**死死依赖点 $k$**（中间那维 $w_k$ 会随 $k$ 变），这就不满足四边形不等式的前提，决策单调性自然也保不住。
唉！所以，$O(n^3)$ 往后**并没有一个便宜的 $O(n^2)$ 可以捡**。
## 正解：Hu-Shing 算法（100pts）
这个算法由 T. C. Hu 与 M. T. Shing 在 1981—1984 年提出，是矩阵链问题目前的最优复杂度做法，并且 $O(n\log n)$ 已被证明是这类问题的下界。它在国内资料很少，但是在韩国竞赛社区被讲得比较透（Baekjoon 上就有 $n\le 2\times 10^5$ 的同题）。顺带一提，它的核心引理（论文里的 Lemma1）最初那版证明其实是**错的**，多年后才由 Shing 本人补上正确证明 —— 一个算法能又快又正确地落地，中间的曲折比想象中多。所以下面我讲清**它在做什么与我的代码是怎么实现的**，严格的最优性证明会在文末给出参考链接，~~因为我太蒻了~~。
### 1.矩阵链就是凸多边形三角剖分
Hu-Shing 的第一记妙招，是把矩阵链问题**转化**成一个几何问题。~~妙！！！！~~
把 $n+1$ 个维度 $w_0,w_1,\dots,w_n$ 看成一个凸多边形的 $n+1$ 个顶点，第 $i$ 个顶点的权值就是 $w_i$。那么：
- 多边形的每条边，对应链里的一个矩阵，而底边对应最终结果。
- 把这个 $(n+1)$ 顶点的多边形用不相交的对角线剖成 $n-1$ 个三角形，**每一种三角剖分方案恰好对应一种加括号方案**。
- 一个顶点为 $w_i,w_j,w_k$ 的三角形，代价就是 $w_i\cdot w_j\cdot w_k$，对应「**把两段子链在中间那个公共维度上乘起来**」。
于是，**最优加括号=让所有三角形权值之和最小的三角剖分**。整个问题被搬到了几何上，接下来就能用~~几何直觉~~去啃它。
### 2.Hu-Shing 三部曲
**（1）找最小顶点。** 找到权值最小的顶点，把它转到 $v_0$ 的位置。直觉上，最便宜的那个维度像一个「**万能因子**」，它出现在越多三角形里越划算，所以以它为顶点拉一把「**扇形**」是很强的候选方案。代码里就是先找最小值下标 $p$，再旋转，最后把多边形闭合。
**（2）扇形分解（利用单调栈一次扫描）。** 一个「**扇形**」指一堆共用同一个顶点（局部最小值）当扇柄的三角形。Hu-Shing 证明了：最优解具有**嵌套扇形**的结构。代码用一个单调栈一次扫过所有顶点，每遇到一个不大于栈顶的顶点，就切出一个扇形单元，记下它的左右边界，并把落在它区间内的、更小的扇形收作「**儿子**」。这些父子关系，就是为后面奠基。
**（3）比值判定 + 两次归约**。一个扇形直接按「**锚在较小端顶点**」来剖不一定最优，有时候应该把某个儿子扇形并进父亲（或反过来）。判断依据是一个**比值** $\dfrac{kp}{kq}$：
- $kp$ 是这一块**分子**，$kq$ 是**分母**，如何求出这个值在代码中会体现。
- **第一次归约**：当**最大**的儿子比值 $\ge \min(v_{l},v_{r})$（较小的那个边界权值）时，把这个儿子吸收进来。
- 然后算出「**本扇形锚在较小端顶点** $g$」时自己的最优代价 $kp=num\cdot v_g$。
- **第二次归约**：当本扇形的比值 $\le$ 剩下最大儿子的比值时，继续把儿子并上来。
这里 $\text{span}(l,r)=\Big(\sum_{k=l+1}^{r} v_k v_{k-1}\Big)-v_l v_r$，这个是**边界折线的代价减去弦**，用前缀和就能轻松 $O(1)$ 查询。
每个扇形要反复取出「**比值最大的儿子**」并合并子结构，这天然是一个**可并堆**能干的活 —— 配对堆，来吧。
### 3.代码里的数据结构必须手写
这份代码之所以长、之所以满屏单字母数组，全是为了 $2\times 10^6$ 这个规模。三处设计，每一处都是被数据逼出来的：
- 利用**侵入式双向链表**存每个扇形的儿子列表：归约时要频繁「**把某个儿子从中间删掉**」，链表 $O(1)$ 删除，且**不额外开内存**。
- 利用**手写配对堆**维护儿子比值：这里必须用**可并堆**。配对堆的 meld 是均摊 $O(1)$，正是我们需要的。
- **前缀和**：$span$ 靠前缀和 $O(1)$；而 $num\cdot v_g$、以及累加的 $ans$ 都可能瞬间超过 long long，需要使用 `__int128`，最后再回到 long long 输出。
还有一个容易被忽略但很致命的点：**空间**。这份代码开了 $13$ 个长度 $2.1\times 10^6$ 的 long long 数组。如果像很多教科书实现那样「**给每个堆结点 $new$ 一下**」，指针开销加上分配器碎片，$2\times 10^6$ 个结点会直接 MLE。**所以还得是手写侵入式结构。**
### 完整代码
[code](https://blog.csdn.net/Tsntsn123/article/details/163528129?spm=1001.2014.3001.5502)
### 复杂度
单调栈一次扫描是 $O(n)$；每个扇形单元至多被弹出、合并常数次，配对堆操作均摊 $O(\log n)$，总时间 $O(n\log n)$，空间 $O(n)$。
### 只想要 $60$ 分？Hu-Shing 还有弱化版
那 $n\le 2\times 10^5$ 的 $60\%$，本质上是**逃不掉 Hu-Shing**的本质思路的。但好消息是在这个规模下，你可以**不用手写**链表和配对堆了！一下这两处可以这样改：
- 儿子列表用 list 或干脆 vector 存下标；
- 比值集合用 priority_queue/multiset，合并时**启发式合并（小的并进大的）**，但是复杂度退化到 $O(n\log^2 n)$。
$n=2\times 10^5$ 时加上 STL 常数也稳过。代价是常数和内存都会有点大 —— 到 $2\times 10^6$ 就会 TLE / MLE，但是核心逻辑与正解**一模一样**，这里就不再贴一份几乎重复的长代码占版面了，把上面那份换成 STL，把小到大合并加上，就是一份稳拿 $60$ 分的实现。
## 结语
这题的价值，其实不在最后那份代码本身，而在这条路：当你写下暴力 $O(n^3)$ 时，你理解了问题，四边形不等式的失败后，你看清了「**转移代价依赖决策点**」这个本质障碍，而 Hu-Shing 用「**几何问题+扇形结构+比值计算**」漂亮地绕了过去。如果这篇题解帮你省下了在题解区枯坐的那几个夜晚，那它就没白写。
> The purpose of computing is insight, not numbers. —— Richard Hamming
祝AC。
## 参考资料
- [洛谷 P1753 矩阵链排序问题](https://www.luogu.com.cn/problem/P1753)。
- [HandWiki, Matrix chain multiplication](https://handwiki.org/wiki/Matrix_chain_multiplication)。
- [node-weighted 多边形三角剖分归约与修正证明的整理](https://arxiv.org/pdf/2104.01777)。
- [On the Parenthesisations of Matrix Chains](https://arxiv.org/pdf/2303.17352)。
- [GitHub, junodeveloper/Hu-Shing](https://github.com/junodeveloper/Hu-Shing)。
- [知乎，矩阵链乘的 O(nlogn) 算法是怎么实现的](https://www.zhihu.com/question/644452577)？
- [EagleBear2002 的博客 P1753 笔记](https://eaglebear2002.github.io/51406/)。
- 《算法导论》。