矩阵优化 DP 以及全局平衡二叉树实现的动态 DP 学习笔记

矩阵乘法

求斐波那契数列的第 \(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} \]

\(f_1\) 和 \(f_2\) 是已知的，后面乘的矩阵是相同的，所以可以用矩阵快速幂。

常见优化矩阵乘法的方法：

预处理出转移矩阵的 \(2\) 的若干次幂次方。

如果转移是 \(1\times n\) 初始矩阵乘以 \(n\times n\) 的转移矩阵的若干次方，可以用初始矩阵依次乘以转移矩阵的 \(2\) 的若干次幂次方，这样一次乘法就由 \(O(n^3)\) 变成 \(O(n^2)\)。

如果转移矩阵很多位置是空，可以把不是零的位置拎出来做乘法。

单位矩阵：只有从左上到右下的对角线的位置是 \(1\) 其它都是 \(0\) 的矩阵。

可以动态修改的矩阵优化 DP

Ex - 01? Queries (atcoder.jp)

[ABC246Ex] 01? Queries - 洛谷 | 计算机科学教育新生态 (luogu.com.cn)

给定长度为 \(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)

P4751 【模板】"动态DP"&动态树分治（加强版） - 洛谷 | 计算机科学教育新生态 (luogu.com.cn)

给你一个 \(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\) 级别。

修改的时候，跳重边可以直接更改父亲的转移矩阵乘积，跳轻边可以更改父亲的转移矩阵。

附上全局平衡二叉树实现上面例题的代码：

点击开 D

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)

容易写出转移方程：\(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} \]

所以只需要记录矩阵的左上角和左下角两个量即可。

点击开 D

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

[P3781 [SDOI2017] 切树游戏 - 洛谷 | 计算机科学教育新生态 (luogu.com.cn)](