动态 DP 学习笔记
不得不承认，去年提高组 D2T3 对动态 DP 起到了良好的普及效果。
动态 DP 主要用于解决一类问题。这类问题一般原本都是较为简单的树上 DP 问题，但是被套上了丧心病狂的修改点权的操作。举个例子，我们来看一道例题。
给定一棵 \(n\) 个点的树。\(i\) 号点的点权为 \(a_i\)。有 \(m\) 次操作，每次操作给定 \(u, w\)，表示修改点 \(u\) 的权值为 \(w\)。你需要在每次操作之后求出这棵树的最大权独立集的权值大小。
我们首先考虑没有修改的情况下怎么做。首先先选取 \(1\) 号点作为全树的根。然后我们设 \(f_{i, 0}\) 表示不选择 \(i\) 号点时，以 \(i\) 号点为根的子树的最大权独立集；\(f_{i, 1}\) 表示选择 \(i\) 号点时，以 \(i\) 号点为根的子树的最大权独立集。我们可以很容易地写出如下的方程：
这里 \(j\) 表示 \(i\) 号点的所有儿子。特殊地，若点 \(i\) 为叶子节点，\(f_{i, 0} = 0, f_{i, 1} = a_i\)。
最后的答案就是 \(\max(f_{1, 0}, f_{1, 1})\)。
接下来带上修改。
首先根据动态规划的转移方程可以发现，我们修改了一个点的点权，只会更改从这个点到根这条路径上节点的 DP 值，其他值是不会发生更改的。这时候如果我们要对整棵树重新求一遍最大权独立集，未免太过浪费。所以我们希望能够更改这条链上的 DP 值。
由于树可能会退化成一条链，这样每次更新就是 \(\mathcal{O(n)}\) 的，显然不可接受。我们希望这条链只更新 \(\log n\) 次……
点分治！抱歉博主太弱了，不会那个被称作“全局平衡二叉树”的厉害做法。
这时候我们请出解决树上问题的神器——重链剖分。
重链剖分有一些性质，这些性质正是它在动态 DP 中能够发挥作用的重要保障。
- 每个点到根的路径上，最多经过 \(\log n\) 条轻边。也就是说，重链的条数最多也只有 \(\log n\) 条。这为动态 DP 的时间复杂度做了保障。
- 每条重链的链尾都是叶子节点，且只有叶子节点没有重儿子。这为动态规划的初始状态和转移方式做了保障。
- 重链剖分中，一条重链所在的区间在剖出的 DFS 序上，是连续的一段区间。这为可以使用数据结构维护区间信息，达到快速转移做了保障。
那么在宏观上，我们相当于在更新时，对于这些重链暴力地互相转移更新。接下来我们考虑一些微观问题：在一条链里，怎么支持快速修改和查询这条链的 DP 值。
我们保持 \(f\) 数组的定义不变。为了迎合重链剖分划分出了轻重儿子，我们形式化地定义 \(g\) 数组：\(g_{i, 1}\) 表示 \(i\) 号点的所有轻儿子，都不取的最大权独立集；\(g_{i, 0}\) 表示 \(i\) 号点的所有轻儿子，可取可不取形成的最大权独立集。这样就可以把上述的 DP 式子大大简化了（至少没有了那个 \(\Sigma\)）。
这里的 \(j\) 表示 \(i\) 号点的重儿子。特殊地，对于叶子节点，\(g_{i, 0} = g_{i, 1} = 0\)。
但是感觉这玩意儿好像不大优美？第二个转移式子中，\(g_{i, 1}\) 和 \(a_i\) 都只和 \(i\) 有关，那么我们不妨把它们合并起来。我们重新定义 \(g_{i, 1}\)：表示 \(i\) 号点只考虑轻儿子的取自己的最大权独立集。那么这时候，第二个方程就可以变为 \(f_{i, 1} = g_{i, 1} + f_{j, 0}\)。
但是这玩意儿咋区间维护嘞？回想一下当初学习斐波那契的时候，我们碰到过这样的 DP 方程：
这个方程涉及上一步的贡献，没法满足结合率，不太舒服。于是我们定义了一个矩阵，化加为乘，于是我们愉快地用快速幂 AC 了。
这道题我们也给它套个矩阵。对于每个点，都表示一个状态，这个状态共有两个值，于是我们考虑维护一个 \(1 \times 2\) 的矩阵。
现在我们要从一个点的重儿子 \(j\) 转移到 \(i\) 上，也就是说我们需要构造出一个转移矩阵使得 \(\begin{bmatrix} f_{j, 0} & f_{j, 1} \end{bmatrix}\) 能够转移到 \(\begin{bmatrix} f_{i, 0} & f_{i, 1} \end{bmatrix}\)。但是我们回顾一下这个转移方程（已更改 \(g_{i, 1}\) 的定义）：
它一点也不满足矩阵乘法的形式啊！
别慌……我们大胆地重定义矩阵乘法！
我们定义一个新的运算符 \(*\)，对于矩阵 \(\mathrm{A}, \mathrm{B}\)，定义 \(\mathrm{A} * \mathrm{B}\) 的结果 \(\mathrm{C}\)，满足：
实现到代码上，就是
struct Matrix {
int mat[MaxN][MaxN];
};
inline Matrix operator * (Matrix a, Matrix b) {
Matrix c;
for (int i = 0; i < n; ++i)
for (int j = 0; j < n; ++j)
for (int k = 0; k < n; ++k)
c.mat[i][j] = max(c.mat[i][j], a.mat[i][k] + b.mat[k][j]);
return c;
}
但是这个东西为什么具有结合率呢？
- 一种感性的理解：由于 \(\max\) 操作和加法操作都是满足结合率的，所以这个运算满足结合率。
- 一种理性但不太严谨的证明：读者不妨拿出之笔，计算几组 \((\mathrm{A} * \mathrm{B}) * \mathrm{C}\) 和 \(\mathrm{A} * (\mathrm{B} * \mathrm{C})\) 的值（如果您计算比较厉害，带上参数算当然更好）。一般情况下，证明了三个满足条件，对于所有情况都是能满足条件的。
于是我们口胡完了结合率的证明。那么我们就可以用了。接下来我们要构造一个转移矩阵，这个是相对难的一个内容。我就介绍一下我个人构造转移矩阵的拙劣方法吧。
在构造一个转移矩阵之前，我们先想办法把这玩意儿变形，变得和运算 \(*\) 差不多。
接着我们把已知的状态和要转移到的状态写在一起，把未知的转移矩阵用 \(\mathrm{U}\) 表示。
我们原来是一个 \(1 \times 2\) 的矩阵，要形成一个 \(1 \times 2\) 的矩阵，那么 \(\mathrm{U}\) 应当是一个 \(2 \times 2\) 的矩阵。那么我们设矩阵左上、右上、左下、右下四个位置分别为 \(u_1, u_2, u_3, u_4\)。接下来把每个位置对应上去。
\(f_{i, 0}\) 的值应该为 \(\max(f_{j, 0} + u_1, f_{j, 1} + u_3)\)。对应转移方程，我们发现 \(u_1\) 应该就是 \(g_{i, 0}\)，\(u_3\) 也是 \(g_{i, 0}\)。同样的，\(f_{i, 1}\) 的值应该为 \(\max(f_{j, 0} + u_2, f_{j, 1} + u_4)\)。对应转移方程，我们发现 \(u_2\) 应该是 \(g_{i, 1}\)，而不存在 \(f_{j, 1}\) 项，就将 \(u_4\) 赋为 \(-\infty\)。最后写出来，检查一遍：
嗯……好像没问题？
这样子，我们对于一条重链，我们的叶子节点就存储了最初始的值，链上每个节点都对应着一个转移矩阵。我们发现这个转移矩阵和重链信息是没有任何关系的，且因为这个矩阵满足结合率，对于一条重链，我们可以之间线段树维护区间乘积（或者叫……“\(*\) 积”？）。然后到了一条重链链头，因为这个点是它父亲的轻儿子，我们需要更新它父亲节点所在的点的转移矩阵。这样子一直跳到根节点就可以了。貌似……大功告成？
重链剖分剖出的 DFS 序，由于先访问了链头，所以这个区间中，链头在区间左端，链尾在区间右端。我们存储的初始信息在叶子节点（也就是链尾）上，因此我们的矩阵 \(*\) 法应当是转移矩阵在前，要维护的值矩阵在后。我们要把这个矩阵前后换个顺序，再转个个儿，加上一些推算，可以变形成：
这样就真的做完了。最后我写一些关于代码实现的小细节：
- 对于一个点查其 dp 值，需要从这个点一直查到区间链尾。因此，树剖时我们需要多维护一个 \(\texttt{End[i]}\)（这里的 \(i\) 是一条重链的链头），表示以 \(i\) 为链头的这条链，链尾（叶子）节点在 DFS 序上的位置。
- 更新线段树上某个点的转移矩阵时，传入的如果是矩阵，递归下去常数太大。一个解决方法是，在线段树外，维护一个矩阵组 \(\texttt{Val[i]}\)，表示每个节点对应的转移矩阵。这样在线段树更新找到对应位置时，直接赋值进来即可。
最后贴上代码。
解释一下变量名：
\(\texttt{Id[i]}\) 表示 \(i\) 号点在 DFS 序中的位置，\(\texttt{Dfn[i]}\) 表示在 DFS 序中下标 \(i\) 的位置对应的是什么点（与 \(\texttt{Id[i]}\) 相反），\(\texttt{Fa[i]}\) 是父亲节点，\(\texttt{Siz[i]}\) 是子树大小，\(\texttt{Dep[i]}\) 是该节点深度（好像没什么用），\(\texttt{Wson[i]}\) 是 \(i\) 号节点的重儿子，\(\texttt{Top[i]}\) 表示 \(i\) 号点所在重链链顶编号。
#include <iostream>
#include <cstdio>
#include <cstring>
using namespace std;
const int MaxN = 100000 + 5, MaxM = 200000 + 5;
const int MaxV = 400000 + 5;
const int INF = 0x7F7F7F7F;
struct Matrix {
int mat[2][2];
Matrix() {
memset(mat, -0x3F, sizeof mat);
}
inline Matrix operator * (Matrix b) {
Matrix c;
for (int i = 0; i < 2; ++i)
for (int j = 0; j < 2; ++j)
for (int k = 0; k < 2; ++k)
c.mat[i][j] = max(c.mat[i][j], mat[i][k] + b.mat[k][j]);
return c;
}
};
int N, M; int cntv, cnte;
int A[MaxN];
int Fa[MaxN], Siz[MaxN], Dep[MaxN], Wson[MaxN];
int Top[MaxN], Id[MaxN], Dfn[MaxN], End[MaxN];
int F[MaxN][2];
int Head[MaxN], To[MaxM], Next[MaxM];
Matrix Val[MaxN];
struct SegTree {
int L[MaxV], R[MaxV];
Matrix M[MaxV];
inline void Push_up(int i) {
M[i] = M[i << 1] * M[i << 1 | 1];
}
void Build_Tree(int left, int right, int i) {
L[i] = left, R[i] = right;
if (L[i] == R[i]) {
M[i] = Val[Dfn[L[i]]];
return;
}
int mid = (L[i] + R[i]) >> 1;
Build_Tree(L[i], mid, i << 1);
Build_Tree(mid + 1, R[i], i << 1 | 1);
Push_up(i);
}
void Update_Tree(int x, int i) {
if (L[i] == R[i]) {
// 直接赋值，减小常数
M[i] = Val[Dfn[x]];
return;
}
int mid = (L[i] + R[i]) >> 1;
if (x <= mid) Update_Tree(x, i << 1);
else Update_Tree(x, i << 1 | 1);
Push_up(i);
}
// 查询一个点的 DP 值，相当于查询这条重链上链尾矩阵和链中转移矩阵的 '*' 积
Matrix Query_Tree(int left, int right, int i) {
if (L[i] == left && R[i] == right) return M[i];
int mid = (L[i] + R[i]) >> 1;
if (right <= mid)
return Query_Tree(left, right, i << 1);
else if (left > mid)
return Query_Tree(left, right, i << 1 | 1);
else
return Query_Tree(left, mid, i << 1) * Query_Tree(mid + 1, right, i << 1 | 1);
}
} T;
inline void add_edge(int from, int to) {
cnte++; To[cnte] = to;
Next[cnte] = Head[from]; Head[from] = cnte;
}
void readin() {
scanf("%d %d", &N, &M);
for (int i = 1; i <= N; ++i)
scanf("%d", &A[i]);
for (int i = 1; i < N; ++i) {
int u, v;
scanf("%d %d", &u, &v);
add_edge(u, v); add_edge(v, u);
}
}
void dfs1(int u) {
Siz[u] = 1;
for (int i = Head[u]; i; i = Next[i]) {
int v = To[i];
if (v == Fa[u]) continue;
Fa[v] = u; Dep[v] = Dep[u] + 1;
dfs1(v);
Siz[u] += Siz[v];
if (Siz[v] > Siz[Wson[u]]) Wson[u] = v;
}
}
void dfs2(int u, int chain) {
cntv++;
Id[u] = cntv; Dfn[cntv] = u;
Top[u] = chain;
End[chain] = max(End[chain], cntv);
// 第二次树剖时直接更新 F, G 数组（这里直接将 G 放入矩阵更新）
F[u][0] = 0, F[u][1] = A[u];
Val[u].mat[0][0] = Val[u].mat[0][1] = 0;
Val[u].mat[1][0] = A[u];
if (Wson[u] != 0) {
dfs2(Wson[u], chain);
// 依照定义，重儿子不应计入 G 数组
F[u][0] += max(F[Wson[u]][0], F[Wson[u]][1]);
F[u][1] += F[Wson[u]][0];
}
for (int i = Head[u]; i; i = Next[i]) {
int v = To[i];
if (v == Fa[u] || v == Wson[u]) continue;
dfs2(v, v);
F[u][0] += max(F[v][0], F[v][1]);
F[u][1] += F[v][0];
Val[u].mat[0][0] += max(F[v][0], F[v][1]);
Val[u].mat[0][1] = Val[u].mat[0][0];
Val[u].mat[1][0] += F[v][0];
}
}
void init() {
readin();
dfs1(1); dfs2(1, 1);
}
void update_path(int u, int w) {
Val[u].mat[1][0] += w - A[u];
A[u] = w;
Matrix bef, aft;
while (u != 0) {
// 计算贡献时，应当用一个 bef 矩阵还原出少掉这个轻儿子的情况，再将 aft 加入更新
bef = T.Query_Tree(Id[Top[u]], End[Top[u]], 1);
T.Update_Tree(Id[u], 1);
aft = T.Query_Tree(Id[Top[u]], End[Top[u]], 1);
u = Fa[Top[u]];
Val[u].mat[0][0] += max(aft.mat[0][0], aft.mat[1][0]) - max(bef.mat[0][0], bef.mat[1][0]);
Val[u].mat[0][1] = Val[u].mat[0][0];
Val[u].mat[1][0] += aft.mat[0][0] - bef.mat[0][0];
}
}
void solve() {
T.Build_Tree(1, N, 1);
for (int i = 1; i <= M; ++i) {
int u, w;
scanf("%d %d", &u, &w);
update_path(u, w);
Matrix Ans = T.Query_Tree(Id[1], End[1], 1);
printf("%d\n", max(Ans.mat[0][0], Ans.mat[1][0]));
}
}
int main() {
init();
solve();
return 0;
}
我们再看一道例题，加深一下对动态 DP 的理解吧。
给定一棵 \(n\) 个点的树，每个点有一个正数点权 \(a_i\)，\(1\) 号节点为根。封住一个点的花费为这个点的点权。要求支持两种操作：
- 将某个点 \(u\) 的点权加上一个正整数 \(w\)；
- 查询封死以 \(x\) 为根的子树的最小花费。封死定义为，对于这棵子树内的所有叶子节点，到 \(x\) 号点的路径上至少要有一个点被封住。
由于 \(w\) 是正数，所以点权只会变大，这个性质使得这题有不用动态 DP 的做法。但是我们仍旧介绍动态 DP 做法。
首先仍然考虑没有修改操作怎么做。设 \(f_i\) 表示封死以 \(i\) 号点为根的子树的最小代价。特殊地，对于叶子节点，\(f_i = a_i\)。
考虑转移。这棵子树，要么在 \(i\) 号点被封，要么它所有儿子子树都被封。我们设 \(j\) 为 \(i\) 的儿子节点，那么有
接下来考虑修改。同样地，我们形式化地定义 \(g_i\) 表示 \(i\) 号节点所有轻儿子的子树都被封上的最小花费。那么如果设 \(j\) 为 \(i\) 的重儿子，该方程可以改写为
注意到这里转移的时候有两种运算，分别是 \(\min\) 和加法运算。这两种运算都是满足结合率的，于是我们定义矩阵新运算 \(*\)：如果 \(\mathrm{A} * \mathrm{B} = \mathrm{C}\)，那么有：
接下来我们考虑如何设置状态矩阵和转移矩阵。
我们注意到转移方程当中的 \(f_j + g_i\) 是一个二项式，这意味着这个矩阵应当是二维的矩阵。可是我们的状态只有一维，这该怎么办呢？
回到这个转移方程上来。我们现在所希望的是，将 \(f_j\) 这个状态转移到 \(f_i\) 状态。也就是说，转移矩阵的值应当只与 \(i\) 有关。那么 \(a_i\) 和 \(g_i\) 这两个值都应该出现在这个转移矩阵上。我们发现 \(g_i\) 已经与 \(f_j\) 对应，而 \(a_i\) 这一项可以看作 \(a_i + 0\)，也就是说状态矩阵另外一维应当是 \(0\)。我们大致可以写出如下的转移方式（未知的地方用 \(u_1, u_2\) 表示）：
为了得出另一个状态矩阵的 \(0\)，对应过来也很容易得出 \(u_1 = \infty, u_2 = 0\)。
这样对于每一个节点，保存一个转移矩阵；特殊地，对于叶子节点，保存状态矩阵，即 \(\begin{bmatrix} a_i \\ 0 \end{bmatrix}\)。然后线段树维护区间信息即可。
同样贴上我丑陋的代码。
#include <iostream>
#include <cstdio>
#include <cstring>
using namespace std;
typedef long long llt;
const int MaxN = 200000 + 5, MaxM = 400000 + 5;
const int MaxV = 800000 + 5;
const llt INF = 0x7F7F7F7F7F7F;
struct Matrix {
llt mat[2][2];
Matrix() {
mat[0][0] = mat[0][1] = mat[1][0] = mat[1][1] = INF;
}
inline Matrix operator * (Matrix b) {
Matrix c;
for (int i = 0; i < 2; ++i)
for (int j = 0; j < 2; ++j)
for (int k = 0; k < 2; ++k)
c.mat[i][j] = min(c.mat[i][j], mat[i][k] + b.mat[k][j]);
return c;
}
};
int N, Q; int cntv, cnte;
llt A[MaxN];
int Head[MaxN], To[MaxM], Next[MaxM];
int Fa[MaxN], Dep[MaxN], Siz[MaxN], Wson[MaxN];
int Dfn[MaxN], Id[MaxN], Top[MaxN], End[MaxN];
llt F[MaxN], G[MaxN];
Matrix Val[MaxN];
struct SegTree {
int L[MaxV], R[MaxV];
Matrix Mul[MaxV];
inline void Push_up(int i) {
Mul[i] = Mul[i << 1] * Mul[i << 1 | 1];
}
void Build_Tree(int left, int right, int i) {
L[i] = left, R[i] = right;
if (L[i] == R[i]) {
Mul[i] = Val[Dfn[left]];
return;
}
int mid = (L[i] + R[i]) >> 1;
Build_Tree(left, mid, i << 1);
Build_Tree(mid + 1, right, i << 1 | 1);
Push_up(i);
}
void Update_Tree(int x, int i) {
if (L[i] == R[i]) {
Mul[i] = Val[Dfn[x]];
return;
}
int mid = (L[i] + R[i]) >> 1;
if (x <= mid) Update_Tree(x, i << 1);
else Update_Tree(x, i << 1 | 1);
Push_up(i);
}
Matrix Query_Tree(int left, int right, int i) {
if (L[i] == left && R[i] == right) return Mul[i];
int mid = (L[i] + R[i]) >> 1;
if (right <= mid) return Query_Tree(left, right, i << 1);
else if (left > mid) return Query_Tree(left, right, i << 1 | 1);
else return Query_Tree(left, mid, i << 1) * Query_Tree(mid + 1, right, i << 1 | 1);
}
} T;
inline char ge_ch() {
char c;
do c = getchar(); while (c != 'C' && c != 'Q');
return c;
}
inline void add_edge(int from, int to) {
cnte++; To[cnte] = to;
Next[cnte] = Head[from]; Head[from] = cnte;
}
void init() {
scanf("%d", &N);
for (int i = 1; i <= N; ++i) scanf("%lld", &A[i]);
for (int i = 1; i < N; ++i) {
int u, v;
scanf("%d %d", &u, &v);
add_edge(u, v); add_edge(v, u);
}
scanf("%d", &Q);
}
void dfs1(int u) {
Siz[u] = 1;
for (int i = Head[u]; i; i = Next[i]) {
int v = To[i];
if (v == Fa[u]) continue;
Fa[v] = u; Dep[v] = Dep[u] + 1;
dfs1(v);
Siz[u] += Siz[v];
if (Siz[v] > Siz[Wson[u]]) Wson[u] = v;
}
}
void dfs2(int u, int chain) {
cntv++;
Dfn[cntv] = u; Id[u] = cntv;
Top[u] = chain; End[chain] = max(End[chain], cntv);
llt sum = 0; F[u] = A[u];
if (Wson[u] != 0) {
dfs2(Wson[u], chain);
sum += F[Wson[u]];
}
for (int i = Head[u]; i; i = Next[i]) {
int v = To[i];
if (v == Fa[u] || v == Wson[u]) continue;
dfs2(v, v);
sum += F[v]; G[u] += F[v];
}
if (Wson[u] != 0) F[u] = min(F[u], sum);
}
void update_path(int u, llt w) {
if (Wson[u] == 0) Val[u].mat[0][0] = w;
else Val[u].mat[0][1] = w;
A[u] = w;
Matrix bef, aft;
while (u != 0) {
bef = T.Query_Tree(Id[Top[u]], End[Top[u]], 1);
T.Update_Tree(Id[u], 1);
aft = T.Query_Tree(Id[Top[u]], End[Top[u]], 1);
u = Fa[Top[u]];
Val[u].mat[0][0] += aft.mat[0][0] - bef.mat[0][0];
}
}
void solve() {
dfs1(1); dfs2(1, 1);
for (int i = 1; i <= N; ++i) {
if (Wson[i] == 0) { // 特殊处理叶子节点
Val[i].mat[0][0] = F[i];
Val[i].mat[1][0] = 0;
continue;
} else {
Val[i].mat[0][0] = G[i];
Val[i].mat[0][1] = A[i];
Val[i].mat[1][0] = INF;
Val[i].mat[1][1] = 0;
}
}
T.Build_Tree(1, N, 1);
for (int i = 1; i <= Q; ++i) {
char opt = ge_ch();
if (opt == 'C') {
int x; llt w;
scanf("%d %lld", &x, &w); w += A[x];
update_path(x, w);
} else {
int x;
scanf("%d", &x);
Matrix Ans = T.Query_Tree(Id[x], End[Top[x]], 1);
printf("%lld\n", Ans.mat[0][0]);
}
}
}
int main() {
init();
solve();
return 0;
}
最后小小的总结一下吧。一个问题一旦支持修改，且如果没有修改的情况下是一个简单的树上 DP 问题，如果实在想不出更简单的做法，可以往动态 DP 的方向去思考。因为动态 DP 的码量巨大，也不大容易调试。以上这几题的代码，我都因为各种各样的小错误而调了很久。在 D2T3 的实战中，我的几位学长们看出来了这是动态 DP，但是迫于时间问题，只写了个暴力上去。所以我个人认为，巩固一些基础的算法、提升自己思维能力这两点，比学习一些高级算法更为重要。保证自己基础暴力分拿到，再去花时间想满分做法，才应当是竞赛中较为优秀的一种做题策略。
完结，撒花~
