矩阵优化 DP 以及全局平衡二叉树实现的动态 DP 学习笔记
fydj
·
·
算法·理论
矩阵乘法
求斐波那契数列的第 n 项，其中 n\le 10^{18}，对数 m 取模。
写出转移方程，f_i=f_{i-1}+f_{i-2}，直接算是 O(n) 的，似乎没什么办法优化。
定义大小分别为 n\times p 和 p\times m 的矩阵（分别为 a 和 b）相乘的结果（矩阵 c）是一个大小为 n\times m 的矩阵，其中 c_{i,j}=\sum_{k=1}^p a_{i,k}\times b_{k,j}。
把转移方程写成矩阵形式，得到：
\begin{bmatrix}
f_i&f_{i-1}
\end{bmatrix}
=
\begin{bmatrix}
f_{i-1}&f_{i-2}
\end{bmatrix}
\times
\begin{bmatrix}
1&1\\1&0
\end{bmatrix}
常见优化矩阵乘法的方法：
预处理出转移矩阵的 $2$ 的若干次幂次方。
如果转移是 $1\times n$ 初始矩阵乘以 $n\times n$ 的转移矩阵的若干次方，可以用初始矩阵依次乘以转移矩阵的 $2$ 的若干次幂次方，这样一次乘法就由 $O(n^3)$ 变成 $O(n^2)$。
如果转移矩阵很多位置是空，可以把不是零的位置拎出来做乘法。
单位矩阵：只有从左上到右下的对角线的位置是 $1$ 其它都是 $0$ 的矩阵。
# 可以动态修改的矩阵优化 DP
[Ex - 01? Queries (atcoder.jp)](https://atcoder.jp/contests/abc246/tasks/abc246_h)
[[ABC246Ex\] 01? Queries - 洛谷 | 计算机科学教育新生态 (luogu.com.cn)](https://www.luogu.com.cn/problem/AT_abc246_h)
给定长度为 $N$ 的仅包含 `0`，`1`，`?` 的字符串 $S$，给定 $Q$ 组询问 $(x_1, c_1), (x_2, c_2), \cdots, (x_q, c_q)$，每次将原字符串中 $x_i$ 位置的字符改为 $c_i$，然后输出 $S$ 有多少种非空子序列，`?` 需任意替换为 `0` 或 `1`。对 $998244353$ 取模。
$1 \le N, Q \le 10^5, 1 \le x_i \le N$。
设 $f_{i,0/1}$ 表示处理到第 $i$ 位，最后一位是 $0$ 还是 $1$ 的方案，得到方程：
$$
f_{i,0}=f_{i-1,0}+[S_i\neq 1](f_{i-1,1}+1)\\
f_{i,1}=f_{i-1,1}+[S_i\neq 0](f_{i-1,0}+1)
$$
注意到有常数项 $1$，不是很好维护，于是再转移矩阵里加多一个常数项 $1$。
$$
\begin{bmatrix}
f_{i,0}&f_{i,1}&1
\end{bmatrix}
=
\begin{bmatrix}
f_{i-1,0}&f_{i-1,1}&1
\end{bmatrix}
\times
\begin{bmatrix}
1&[S_i\neq0]&0\\
[S_i\neq1]&1&0\\
[S_i\neq1]&[S_i\neq0]&1
\end{bmatrix}
$$
要求动态修改一个位置的值，相当于修改一个位置的转移矩阵，可以用线段树维护。
时间复杂度 $O(27(N+Q)\log N)$，前面的 $27$ 是矩阵乘法的常数。
# 动态 DP
[P4719 【模板】"动态 DP"&动态树分治 - 洛谷 | 计算机科学教育新生态 (luogu.com.cn)](https://www.luogu.com.cn/problem/P4719)
[P4751 【模板】"动态DP"&动态树分治（加强版） - 洛谷 | 计算机科学教育新生态 (luogu.com.cn)](https://www.luogu.com.cn/problem/P4751)
给你一个 $n$ 个点的有点权（点 $i$ 的点权是 $a_i$）的树，$q$ 次操作，每次修正一个点的点权，每次修正后输出最大权独立集的权值大小。
设 $f_{i,0/1}$ 表示考虑 $i$ 的子树，选或不选点 $i$ 形成的最大权独立集的权值大小，列出 DP 式子：
$$
f_{i,0}=\sum _{j\in son_i} \max(f_{j,0},f_{j,1})\\
f_{i,1}=\sum _{j\in son_i} f_{j,0}+a_i
$$
考虑树链剖分，设 $mson_i$ 为点 $i$ 的重儿子，$lson_i$ 为点 $i$ 的轻儿子，又设：
$$
g_{i,0}=\sum _{j\in lson_i} \max(f_{j,0},f_{j,1})\\
g_{i,1}=\sum _{j\in lson_i} f_{j,0}+a_i
$$
那么可以列出转移方程：
$$
f_{i,0}=g_{i,0}+\max(f_{mson_i,0},f_{mson_i,1})\\
f_{i,1}=g_{i,1}+f_{mson_i,1}
$$
重新定义矩阵乘法：定义大小分别为 $n\times p$ 和 $p\times m$ 的矩阵（分别为 $a$ 和 $b$）相乘的结果（矩阵 $c$）是一个大小为 $n\times m$ 的矩阵，其中 $c_{i,j}=\max_{k=1}^p a_{i,k}+b_{k,j}$。就是把 $\sum$ 换成 $\max$，把 $\times$ 换成 $+$，称这种矩阵乘法为 $\max +$ 矩阵乘法。
那么化成矩阵就是：
$$
\begin{bmatrix}
f_{i,0}\\
f_{i,1}
\end{bmatrix}
=
\begin{bmatrix}
g_{i,0}&g_{i,0}\\
-\inf&g_{i,1}
\end{bmatrix}
\times
\begin{bmatrix}
f_{mson_i,0}\\
f_{mson_i,1}
\end{bmatrix}
$$
叶子（LeaF）节点没有儿子，所以后面乘的是 $\begin{bmatrix}0\\0\end{bmatrix}$。
树链剖分，用线段树维护一条重链上转移矩阵的值，跳轻边就把 $f$ 贡献到父亲的 $g$ 里面，时间复杂度是 $O(n\log^2 n)$。
但你会被卡树剖~~，而且树剖难写~~。
# 全局平衡二叉树
LCT 由于有 splay 的良好性质，所以高度是 $\log n$ 级别的，它可以很好地在 $O(\log n)$ 的时间内维护一条重链上转移矩阵的乘积。
但 LCT 常数又大有难写，注意到树的形态是固定的，用 LCT 会不会大材小用？
有一个叫全局平衡二叉树的一个东西，是一个静态的 LCT，给重链上的点附上一个权值（这个权值只被用来构建全局平衡二叉树），权值是当前点的子树大小减去重儿子的子树大小。一条重链的一段缩成一颗 splay 的时候，找到它的带权中点，以它为根。这样构建出的单棵 splay 可能不一定特别平衡，但从整体来看，整棵树是平衡的，深度不会超过 $\log n$ 级别。
修改的时候，跳重边可以直接更改父亲的转移矩阵乘积，跳轻边可以更改父亲的转移矩阵。
附上全局平衡二叉树实现上面例题的代码：
```cpp
const int N=1000099; const int inf=1047483646;
struct ju { int a[2][2]={}; friend ju operator * (const ju a,const ju b) { ju ans; int i,j; for(i=0;i<2;++i) for(j=0;j<2;++j) ans.a[i][j]=max(a.a[i][0]+b.a[0][j],a.a[i][1]+b.a[1][j]); return ans; } } val[N]={},sum[N]={};
int n,q,a[N]={},w[N]={},sz[N]={},lsz[N]={};
int mson[N]={},type[N]={},stype=0; vector<int> lian[N]={};
int fa[N]={},ch[N][2]={},f[N][2]={},g[N][2]={};
int em=0,e[N*2]={},ls[N]={},nx[N*2]={};
void insert(int x,int y) { e[++em]=y,nx[em]=ls[x],ls[x]=em; e[++em]=x,nx[em]=ls[y],ls[y]=em; return ; }
#define lson (ch[x][0])
#define rson (ch[x][1])
#define get(x) (ch[fa[x]][1]==(x))
#define isroot(x) (ch[fa[x]][0]!=(x)&&ch[fa[x]][1]!=(x))
int makebst(int tp,int l,int r) {
int i,tot=0,s=0,x; for(i=l;x=lian[tp][i],i<=r;++i) tot+=lsz[x];
for(i=l;x=lian[tp][i],s+=lsz[x];++i)
if(s*2>=tot) {
if(l<i) fa[lson=makebst(tp,l,i-1)]=x;
if(i<r) fa[rson=makebst(tp,i+1,r)]=x;
sum[x]=sum[lson]*val[x]*sum[rson];
return x;
}
}
int make(int rt) {
int i,tp=type[rt],y,oldf=fa[rt];
for(auto x:lian[tp]) {
for(i=ls[x];y=e[i],i;i=nx[i])
if(y!=fa[x]&&y!=mson[x])
y=make(y),g[x][0]+=max(f[y][0],f[y][1]),g[x][1]+=f[y][0];
g[x][1]+=a[x],val[x]=(ju){g[x][0],g[x][0],g[x][1],-inf};
}
fa[rt=makebst(tp,0,lian[tp].size()-1)]=oldf;
f[rt][0]=max(sum[rt].a[0][0],sum[rt].a[0][1]);
f[rt][1]=max(sum[rt].a[1][0],sum[rt].a[1][1]);
return rt;
}
void change(int x) {
sum[x]=sum[lson]*val[x]*sum[rson];
if(fa[x]&&isroot(x))
g[fa[x]][0]-=max(f[x][0],f[x][1]),g[fa[x]][1]-=f[x][0];
if(isroot(x))
f[x][0]=max(sum[x].a[0][0],sum[x].a[0][1]),f[x][1]=max(sum[x].a[1][0],sum[x].a[1][1]);
if(fa[x]&&isroot(x))
g[fa[x]][0]+=max(f[x][0],f[x][1]),g[fa[x]][1]+=f[x][0],
val[fa[x]]=(ju){g[fa[x]][0],g[fa[x]][0],g[fa[x]][1],-inf};
if(fa[x]) change(fa[x]);
return ;
}
int main()
{
// usefile("ddp");
int i,x,y,l,r,root,lstans=0; ju nxt;
read(n,q); for(i=1;i<=n;++i) read(a[i]); for(i=1;i<n;++i) read(x,y),insert(x,y);
w[l=r=1]=1; while(l<=r) { x=w[l]; for(i=ls[x];y=e[i],i;i=nx[i]) if(y!=fa[x]) fa[y]=x,w[++r]=y; ++l; }
for(i=n;x=w[i],y=fa[x],i;--i) { ++sz[x],++lsz[x],sz[y]+=sz[x],lsz[y]+=sz[x]; if(!mson[y]||sz[mson[y]]<sz[x]) lsz[y]+=sz[mson[y]]-sz[x],mson[y]=x; }
for(i=1;x=w[i],y=mson[x],i<=n;++i) { if(!type[x]) type[x]=++stype,lian[stype].push_back(x); if(y) type[y]=type[x],lian[type[x]].push_back(y); }
val[0]=sum[0]={0,-inf,-inf,0},root=make(1);
loop : --q; read(x,y),x^=lstans,val[x].a[1][0]+=y-a[x],g[x][1]+=y-a[x],a[x]=y,change(x),printf("%d\n",lstans=max(f[root][0],f[root][1])); if(q) goto loop; return 0;
}
```
# 其它技巧
[P6021 洪水 - 洛谷 | 计算机科学教育新生态 (luogu.com.cn)](https://www.luogu.com.cn/problem/P6021)
容易写出转移方程：$f_i=\min(\sum_{j\in son_i} f_j + a_i)$。
设 $s_i=\sum_{j\in lson_i} f_j$，可以列成矩阵形式，这里的 $\times$ 是 $\min+$ 矩阵乘法：
$$
\begin{bmatrix}
f_i&0
\end{bmatrix}
=
\begin{bmatrix}
f_{mson_i}&0
\end{bmatrix}
\times
\begin{bmatrix}
s_i&\inf\\
a_i&0
\end{bmatrix}
$$
用全局平衡二叉树维护。
这里的转移和上一题不一样，上一题的初始矩阵是列向量，是竖着的，用根节点的转移矩阵乘以重儿子的矩阵，从上乘到下；而这里的初始矩阵是行向量，是横着的，用重儿子的矩阵乘以根节点的转移矩阵，从下乘到上。
但，题目要求询问子树，这相当于询问一条重链的后缀乘积，可以在全局平衡二叉树的 splay 中实现。
时间复杂度 $O(n\log n)$。
你还可以继续卡常，注意到：
$$
\begin{bmatrix}
a&\inf\\
b&0
\end{bmatrix}
\times
\begin{bmatrix}
c&\inf\\
d&0
\end{bmatrix}
=
\begin{bmatrix}
a+c&\inf\\
\min(b+c,d)&0
\end{bmatrix}
$$
所以只需要记录矩阵的左上角和左下角两个量即可。
```cpp
const int N=200099;
int n,q,fa[N]={},sz[N]={},mson[N]={},tag[N]={};
int type[N]={},stype=0,num[N]={},leng[N]={},root[N]={};
vector<int> lian[N]={};
int em=0,e[N*2]={},ls[N]={},nx[N*2]={};
ll a[N]={};
void insert(int x,int y) { e[++em]=y,nx[em]=ls[x],ls[x]=em ;e[++em]=x,nx[em]=ls[y],ls[y]=em; return ; }
struct ju {
ll a,b;
ju(ll a_=0,ll b_=0) { a=a_,b=b_; return ; }
friend ju operator * (const ju a,const ju b) {
return ju(a.a+b.a,min(a.b+b.a,b.b));
}
} v[N]={},sum[N]={}; int lson[N]={},rson[N]={};
void setup(int x) {
int i,y; sz[x]=1;
for(i=ls[x];y=e[i],i;i=nx[i])
if(y!=fa[x]) {
fa[y]=x,setup(y);
sz[x]+=sz[y];
if(sz[mson[x]]<sz[y])
mson[x]=y;
}
tag[x]=sz[x]-sz[mson[x]];
return ;
}
void getlian(int x) {
if(!type[x])
type[x]=++stype,
num[x]=0,leng[stype]=1,
lian[type[x]].push_back(x);
if(!mson[x]) return ;
type[mson[x]]=type[x];
num[mson[x]]=leng[type[x]]++;
lian[type[x]].push_back(mson[x]);
getlian(mson[x]);
int i,y;
for(i=ls[x];y=e[i],i;i=nx[i])
if(y!=fa[x]&&y!=mson[x])
getlian(y);
return ;
}
void update(int x) { sum[x]=sum[rson[x]]*v[x]*sum[lson[x]]; return ; }
int makebst(int tp,int l,int r) {
int i,s=0,tot=0;
for(i=l;i<=r;++i) tot+=tag[lian[tp][i]];
for(i=l;i<=r;++i) {
s+=tag[lian[tp][i]];
if(s*2>=tot) {
int x=lian[tp][i]; v[x].b=a[x];
if(l<i) fa[lson[x]=makebst(tp,l,i-1)]=x;
if(i<r) fa[rson[x]=makebst(tp,i+1,r)]=x;
update(x);
return x;
}
}
}
int make(int x) {
int i,y;
for(auto p:lian[type[x]])
for(i=ls[p];y=e[i],i;i=nx[i])
if(y!=fa[p]&&y!=mson[p])
y=make(y),fa[y]=p,v[p].a+=sum[y].b;
return root[type[x]]=makebst(type[x],0,leng[type[x]]-1);
}
int getopt() { static char chart; do chart=getchar(); while(chart!='Q'&&chart!='C'); return chart=='Q'; }
ju query(int x,int p) {
if(!x) return ju(0,1e15);
else if(num[x]==p)
return sum[rson[x]]*v[x];
else if(num[x]<p)
return query(rson[x],p);
else return sum[rson[x]]*v[x]*query(lson[x],p);
}
void change(int x) {
if(lson[fa[x]]==x||rson[fa[x]]==x)
update(x),change(fa[x]);
else if(!fa[x])
update(x);
else
v[fa[x]].a-=sum[x].b,
update(x),
v[fa[x]].a+=sum[x].b,
change(fa[x]);
return ;
}
int main()
{
// usefile("flood");
int i,x,y,opt;
read(n);
for(i=1;i<=n;++i)
read(a[i]);
for(i=1;i<n;++i)
read(x,y),insert(x,y);
setup(1),getlian(1);
v[0]=sum[0]=ju(0,1e15),fa[make(1)]=0;
read(q);
loop : --q;
opt=getopt();
if(opt) read(x),printf("%lld\n",query(root[type[x]],num[x]).b);
else read(x,y),a[x]+=y,v[x].b+=y,change(x);
if(q) goto loop;
return 0;
}
```
# 容易写错的地方
这里讲述的是个人经验，你可能没有容易写错的地方。
+ 注意构建全局平衡二叉树的时候对父亲的更改，有三个地方要注意：全局平衡二叉树内左右儿子的父亲重新设为根节点；splay 的父亲要设为对应的重链最顶上的节点的父亲；整颗全局平衡二叉树的根节点的父亲要设为 0。
+ 更改的时候要分三类：跳重边，直接更改 splay 子树内转移矩阵的乘积，并更新父亲；跳轻边，先把对父亲转移矩阵的贡献删去，再更改 splay 子树内转移矩阵的乘积以求出 DP 值，然后再把对父亲转移矩阵的贡献加回来，并更新父亲；根节点，直接更改乘积，不需要也不可以更新父亲。
+ $0$ 号节点的转移矩阵和转移矩阵乘积要设为单位矩阵，普通矩阵乘法是 $\begin{bmatrix}1&0\\0&1\end{bmatrix}$，$\min +$ 矩阵乘法是 $\begin{bmatrix}0&\inf\\\inf&0\end{bmatrix}$，若是 $\max +$ 矩阵乘法则把 $\inf$ 改为 $-\inf$。
# 练习
[Problem - 573D - Codeforces](https://codeforces.com/problemset/problem/573/D)
[P3781 [SDOI2017\] 切树游戏 - 洛谷 | 计算机科学教育新生态 (luogu.com.cn)](https://www.luogu.com.cn/problem/P3781)