广义矩阵乘法与动态 DP 学习笔记
我也不知道为什么来学这个，挺有趣的。
通过广义矩阵乘法我们可以通过数据结构维护某一类型的 DP。
广义矩阵乘法
一些约定：
- \(\otimes\) 有交换律：\(a\otimes b=b\otimes a\)；
- \(\otimes\) 有结合律：\((a\otimes b)\otimes c=a\otimes (b\otimes c)\)；
- \(\otimes\) 对 \(\oplus\) 有分配律：\(a\otimes (b\oplus c)=(a\otimes b)\oplus(a\otimes c)\)；
若 \(\otimes\) 满足交换律、结合律，\(\otimes\) 对 \(\oplus\) 有分配律，那么我们就有广义矩阵乘法：
广义矩阵乘法也同样满足普通矩阵乘法的结合律。
常见的有：
- \(\oplus\) 是 \(+\)，\(\otimes\) 是 \(\times\)（普通矩阵乘法）；
- \(\oplus\) 是 \(\min\) 或 \(\max\)，\(\otimes\) 是 \(+\)；
应用
对于这样的一类 DP，以“最大子段和”为例：
给出序列 \(\{a_i\}\)，求一个区间使得区间权值和最大。
设 \(f_i\) 表示以 \(i\) 结尾的最大子段和，那么
答案为 \(ans=\max\{f_i\}\)。
现在把它化成一个好看的形式，记 \(ans_i\) 表示 \(f\) 的前缀最大值：
那么可以把它化为一个 \(\oplus\) 是 \(\max\)，\(\otimes\) 是 \(+\) 的矩阵乘法，这里需要多维护一个 \(0\) 方便转移，
这么做有什么优点么？多了 \(27\) 倍常数！那就是 DP 有了结合律，这样就可以满足很多数据结构的要求了。
如果将问题扩展到区间最大子段和，那么只需要用一个线段树维护区间乘积，或者直接倍增，会好写很多。
这样还可以方便的加上修改操作，这就引出了动态 DP。
动态 DP
先看一道模板题 P4719 【模板】"动态 DP"&动态树分治：
给出 \(n\) 点的树，每个点有权值。\(m\) 次询问，单点修改权值并询问最大权独立集的权值和。
\(1\le n,m\le 10^5\)。
回想静态的最大权独立集怎么求，设 \(f(u,1/0)\) 表示\(u\) 选 / 不选的情况下，以 \(u\) 为根的子树的最大权独立集。
这个看起来和“最大子段和”的问题很像，尝试去转化转移形式。
但是 \(u\) 有多个儿子，按照树上问题的套路，只关心重儿子 \(w_u\)，
这里 \(g(u,1/0)\) 代表轻儿子部分的答案，即
好了可以化成矩阵乘法了，
怎么维护这个东西？因为 \(u\) 的答案是只由重儿子 \(w_u\) 转移而来的，考虑用一个重链剖分维护树上链的乘积，最后答案就是根所在链的乘积。
修改的话，就是修改 \(g(u,1/0)\) 构成的那个矩阵，并且只会影响 \(\log n\) 条重链的乘积。复杂度 \(\mathcal O(m\log^2n)\)。
觉得这题写的挺优雅的，贴个代码。
code
const int N = 2e5 + 5;
const int Inf = 0x3f3f3f3f;
struct Matrix {
int a[2][2];
int* operator [](const int k) { return a[k]; }
const int* operator [](const int k) const { return a[k]; }
Matrix() { a[0][0] = a[0][1] = a[1][0] = a[1][1] = -Inf; }
friend Matrix operator *(const Matrix &x, const Matrix &y) {
Matrix res;
res[0][0] = max(x[0][0] + y[0][0], x[0][1] + y[1][0]);
res[0][1] = max(x[0][0] + y[0][1], x[0][1] + y[1][1]);
res[1][0] = max(x[1][0] + y[0][0], x[1][1] + y[1][0]);
res[1][1] = max(x[1][0] + y[0][1], x[1][1] + y[1][1]);
return res;
}
};
struct Graph {
int tot;
int head[N], nxt[N], to[N];
inline void add(int u, int v) {
++tot, nxt[tot] = head[u], to[tot] = v, head[u] = tot;
++tot, nxt[tot] = head[v], to[tot] = u, head[v] = tot;
}
} G;
int n, Q;
int val[N], f[N][2];
Matrix g[N];
int dfc;
int fa[N], top[N], size[N], son[N], dfn[N], end[N];
void dfs1(int u, int fat) {
fa[u] = fat;
size[u] = 1;
for (int i = G.head[u]; i; i = G.nxt[i]) {
int v = G.to[i];
if (v == fat) continue;
dfs1(v, u);
size[u] += size[v];
if (!son[u] || size[v] > size[son[u]]) son[u] = v;
}
}
void dfs2(int u, int tp) {
dfn[u] = ++dfc;
top[u] = tp;
end[tp] = dfn[u];
f[u][1] = val[u];
if (son[u]) {
int v = son[u];
dfs2(v, tp);
f[u][0] += max(f[v][0], f[v][1]);
f[u][1] += f[v][0];
}
int g0 = 0, g1 = val[u];
for (int i = G.head[u]; i; i = G.nxt[i]) {
int v = G.to[i];
if (v == fa[u] || v == son[u]) continue;
dfs2(v, v);
f[u][0] += max(f[v][0], f[v][1]);
f[u][1] += f[v][0];
g0 += max(f[v][0], f[v][1]);
g1 += f[v][0];
}
Matrix &mat = g[dfn[u]];
mat[0][0] = mat[0][1] = g0;
mat[1][0] = g1;
}
struct SegmentTree {
#define ls(_) ((_) << 1)
#define rs(_) ((_) << 1 | 1)
Matrix f[N*4];
void build(int u, int l, int r) {
if (l == r) return f[u] = g[l], void();
int mid = (l + r) >> 1;
build(ls(u), l, mid), build(rs(u), mid+1, r);
f[u] = f[ls(u)] * f[rs(u)];
}
void update(int u, int l, int r, int p) {
if (l == r) return f[u] = g[p], void();
int mid = (l + r) >> 1;
if (p <= mid) update(ls(u), l, mid, p);
else update(rs(u), mid+1, r, p);
f[u] = f[ls(u)] * f[rs(u)];
}
Matrix query(int u, int l, int r, int ql, int qr) {
if (ql <= l && r <= qr) return f[u];
int mid = (l + r) >> 1;
if (ql > mid) return query(rs(u), mid+1, r, ql, qr);
if (qr <= mid) return query(ls(u), l, mid, ql, qr);
return query(ls(u), l, mid, ql, qr) * query(rs(u), mid+1, r, ql, qr);
}
#undef ls
#undef rs
} T;
inline void modify(int u, int k) {
g[dfn[u]][1][0] += k - val[u];
val[u] = k;
while (u) {
int tp = top[u];
Matrix sav = T.query(1, 1, n, dfn[tp], end[tp]);
T.update(1, 1, n, dfn[u]);
Matrix now = T.query(1, 1, n, dfn[tp], end[tp]);
int v = fa[top[u]];
Matrix &mat = g[dfn[v]];
if (v) {
mat[0][0] += max(now[0][0], now[1][0]) - max(sav[0][0], sav[1][0]);
mat[0][1] = mat[0][0];
mat[1][0] += now[0][0] - sav[0][0];
}
u = v;
}
}
inline void main() {
cin >> n >> Q;
for (int i = 1; i <= n; ++i) cin >> val[i];
for (int i = 1; i < n; ++i) {
int u, v;
cin >> u >> v;
G.add(u, v);
}
dfs1(1, 0), dfs2(1, 1), T.build(1, 1, n);
while (Q--) {
int x, y;
cin >> x >> y;
modify(x, y);
Matrix ans = T.query(1, 1, n, dfn[1], end[1]);
cout << max(ans[0][0], ans[1][0]) << '\n';
}
}
全局平衡二叉树
其实还可以通过 LCT \(\mathcal O(m\log n)\) 做，但是常数很大。
注意到这个 LCT 不需要 link、cut，尝试去静态化它的结构。在最开始以最平衡的形态构建 LCT，之后边不再改变，这样就产生了全局平衡二叉树。
具体的构建方法，先把树重链剖分，然后对于每条重链，加权二分重心，然后递归构建。链上每个点的权值是该点加上其轻儿子的大小。而重链所分出的轻链，就像 LCT 一样只保留父亲关系。
这样相当于每条重链维护了一棵 BST，整棵树则由 BST 森林构成，整体上是平衡的。在这棵树上我们每次只能操作从某一点到根的路径，从轻链跳到重链至多有 \(\log n\) 次；在一条链上的 BST 中向上跳，每次 \(size\) 至少翻倍，而 BST 森林的总 \(size\) 为 \(n\)，所以这部分总步数也是至多 \(\log n\) 次。我们得到每次操作复杂度 \(\mathcal O(\log n)\)。
全局平衡二叉树在很多时候可以代替操作简单的重链剖分，复杂度和常数都比重链剖分优秀。但是个人认为比链剖难写，作用最大的地方只有树上的动态 DP 吧。
回到模板题加强版 P4751 【模板】"动态DP"&动态树分治（加强版）。全局平衡二叉树的做法就是在每条链的 BST 维护子树的矩阵乘积，复杂度 \(\mathcal O(m\log n)\)。
这个代码可能有点丑。
code
// 省略 fread / fwrite
using IO::read;
using IO::print;
const int N = 2e6 + 5;
const int Inf = 0x3f3f3f3f;
struct Matrix {
int a[2][2];
int* operator [](const int k) { return a[k]; }
const int* operator [](const int k) const { return a[k]; }
Matrix() { a[0][0] = a[0][1] = a[1][0] = a[1][1] = -Inf; }
friend Matrix operator *(const Matrix &x, const Matrix &y) {
Matrix res;
res[0][0] = max(x[0][0] + y[0][0], x[0][1] + y[1][0]);
res[0][1] = max(x[0][0] + y[0][1], x[0][1] + y[1][1]);
res[1][0] = max(x[1][0] + y[0][0], x[1][1] + y[1][0]);
res[1][1] = max(x[1][0] + y[0][1], x[1][1] + y[1][1]);
return res;
}
};
struct Graph {
int tot;
int head[N], nxt[N], to[N];
inline void add(int u, int v) {
++tot, nxt[tot] = head[u], to[tot] = v, head[u] = tot;
++tot, nxt[tot] = head[v], to[tot] = u, head[v] = tot;
}
} G;
int n, Q;
int val[N], f[N][2];
Matrix g[N];
int size[N], son[N], fat[N];
void dfs1(int u, int pre) {
fat[u] = pre;
size[u] = 1;
for (int i = G.head[u]; i; i = G.nxt[i]) {
int v = G.to[i];
if (v == pre) continue;
dfs1(v, u);
size[u] += size[v];
if (!son[u] || size[v] > size[son[u]]) son[u] = v;
}
}
void dfs2(int u) {
f[u][0] = 0, f[u][1] = val[u];
int g0 = 0, g1 = val[u];
for (int i = G.head[u]; i; i = G.nxt[i]) {
int v = G.to[i];
if (v == fat[u]) continue;
dfs2(v);
f[u][0] += max(f[v][0], f[v][1]);
f[u][1] += f[v][0];
if (v != son[u]) {
g0 += max(f[v][0], f[v][1]);
g1 += f[v][0];
}
}
g[u][0][0] = g[u][0][1] = g0;
g[u][1][0] = g1;
}
namespace GBT {
int root;
int fa[N], ls[N], rs[N];
int lst[N], sz[N];
bool is_root[N];
Matrix s[N];
inline void push_up(int u) {
s[u] = g[u];
if (ls[u]) s[u] = s[ls[u]] * s[u];
if (rs[u]) s[u] = s[u] * s[rs[u]];
}
int build(int pl, int pr, int fat) {
if (pl >= pr) return 0;
int l = pl, r = pr, mid;
while (r - l > 1) {
mid = (l + r) >> 1;
if (2 * (sz[mid] - sz[pl]) < sz[pr] - sz[pl]) l = mid;
else r = mid;
}
int x = lst[r];
fa[x] = fat;
ls[x] = build(pl, r - 1, x);
rs[x] = build(r, pr, x);
push_up(x);
return x;
}
int divide(int u) {
for (int x = u; x; x = son[x]) {
for (int i = G.head[x]; i; i = G.nxt[i]) {
int y = G.to[i];
if (y == fat[x] || y == son[x]) continue;
fa[divide(y)] = x;
}
}
int cnt = 0;
for (int x = u; x; x = son[x]) {
lst[++cnt] = x;
sz[cnt] = sz[cnt - 1] + size[x] - size[son[x]];
}
int x = build(0, cnt, 0);
is_root[x] = true;
return x;
}
inline void modify(int u, int k) {
g[u][1][0] += k - val[u];
val[u] = k;
while (u) {
if (is_root[u] && fa[u]) {
int v = fa[u];
int sav0 = s[u][0][0], sav1 = s[u][1][0];
push_up(u);
int now0 = s[u][0][0], now1 = s[u][1][0];
Matrix &to = g[v];
to[0][0] += max(now0, now1) - max(sav0, sav1);
to[0][1] = to[0][0];
to[1][0] += now0 - sav0;
}
else push_up(u);
u = fa[u];
}
}
inline int get_ans() {
return max(s[root][0][0], s[root][1][0]);
}
}
inline void main() {
read(n), read(Q);
for (int i = 1; i <= n; ++i) read(val[i]);
for (int i = 1; i < n; ++i) {
int u, v;
read(u), read(v);
G.add(u, v);
}
dfs1(1, 0), dfs2(1);
GBT::root = GBT::divide(1);
int last_ans = 0;
while (Q--) {
int x, y;
read(x), read(y);
x ^= last_ans;
GBT::modify(x, y);
last_ans = GBT::get_ans();
print(last_ans);
}
}
一些题目
P5024 [NOIP2018 提高组] 保卫王国
这题 DDP 不知道比倍增好写多少倍。
先不考虑修改，列出一个 DP，\(f(u,1/0)\) 表示 \(u\) 选 / 不选的最小覆盖代价。
这和最大权独立集几乎一样，化成矩阵乘法（这里 \(\oplus\) 是 \(\min\)，\(\otimes\) 是 \(+\)），
加上修改，如果要求 \(x\) 必须选，那就将 \(a_x\) 改为 \(-\infty\)；要求 \(x\) 不能选，那就将 \(a_x\) 改为 \(+\infty\)。然后就把这题水过去了。
具体实现的时候要注意，若矩阵乘积为 \(A\)，那么 \(A_{0,1},A_{1,1}\) 分别表示 \(f(u,0),f(u,1)\)。
P6021 洪水
也是个板子题。
S2OJ#347 越野赛车问题
link 这是校内 OJ，我还是复述一下题意吧。
\(n\) 个点的树，每条边有限制 \(l_i,r_i\)。每次询问给出 \(k\)，在满足 \(l_i\le k\le r_i\) 的所有边构成的子图中，求最长的简单路径。
\(1\le n\le 7\times 10^4\)，\(1\le l_i\le r_i\le n\)。
这题正解是线段树分治 + 可撤销并查集。不过也可以动态 DP。
类似线段树分治的思路，将询问离线，以 \(l_i,r_i\) 为时间轴处理，可以发现这就是一个动态加删边、求最长直径的问题。
树的直径是可以 DP 的，设 \(f(u,0)\) 表示 \(u\) 为一端的最长链，\(f(u,1)\) 为 \(u\) 为根的子树的答案，
改成可以矩阵转移的形式，
这里需要用一个可删堆来维护 \(g(u,0)\), \(g(u,1)\)，就能支持修改。
得到矩阵
就可以做了。
太难码了，我就口胡一下，不过神仙 \(\mathrm{\color{black}{Y}\color{red}{CS\_GG}}\) 写了代码，看 这里。
P3781 [SDOI2017]切树游戏
大毒瘤，看懂写完调完花了将近一天。
\(\newcommand{\xor}{\operatorname{xor}}\)
设集合幂级数 \(f_u\) 为 \(u\) 为根的子树中答案（\(f_u[k]\) 表示深度最小的点为 \(u\) 且权值为 \(k\) 的连通块数），在一开始 \(f_{u*}=x^{val_u}\)，那么转移为
考虑用 FWT 优化其中的 \(\xor\) 卷积，
注意这里的 \(\prod\) 的乘法是按位相乘，所以 \(1\) 其实相当于一个系数全为 \(1\) 的集合集合幂级数，
最后的答案就是对 \(f_u\) 求和，再维护一个 \(s_u\)，
按照套路分为重儿子和轻儿子转移，设
这样可以用矩阵转移，
但是我们发现用于转移的这个矩阵，第三行和第二列都是永远不变的，
所以实现时矩阵只需要记录和处理四个量，常数由 \(27\) 降为 \(4\)。
复杂度 \(\mathcal O(qm\log n)\)。
参考资料
- 「笔记」广义矩阵乘法与 DP - Luckyblock
- 动态DP基础 - GKxx
浙公网安备 33010602011771号