【学习笔记】动态DP

动态DP学习笔记

动态DP是一种支持DP修改的黑科技，尽管它的名字听起来很别扭。（动态的动态规划）

一、引入

首先我们先来回顾一道经典的DP题：最大子段和。

对于一个序列 \(a\) ，对于所有 \(s=\displaystyle\sum_{i=l}^{r}a_i\) ，求 \(\max s\)。

我们设 \(f_i\) 表示以 \(a_i\) 结尾的最大子段和，有

\[f_i=\max\{f_{i-1}+a_{i},a_i\} \]

最终答案是 \(\max f_i\)。

那如果我们要对 \(a_i\) 进行修改呢？当然可以 \(O(nm)\) 地去求解，但这显然不能满足我们对效率的需求。这就需要用到动态DP。

二、广义矩阵乘法

说到这我们便要引入一个新的运算——广义矩阵乘法。

和普通的矩阵乘法规则有所不同，所以就算你没有学过矩阵也没关系。

（1）定义

定义广义矩阵乘法为：

\[C=A \times B \\ C_{i,j}=\displaystyle\max_{k=1}^{n}({A_{i,k}+B_{k,j}}) \]

看不懂？没关系，说人话就是，两个矩阵的广义相乘得到的矩阵，第 \(i\) 行第 \(j\) 列等于矩阵 \(A\) 第 \(i\) 行与矩阵 \(B\) 第 \(j\) 列每项两两相加的最大值。

还是看不懂？那来算一算吧，比如：

\[\begin{bmatrix} 1 &1\\ 1 &2 \end{bmatrix} \times \begin{bmatrix} 3 &0\\ 1 &1 \end{bmatrix}= \begin{bmatrix} 4 &2\\ 4 &3 \end{bmatrix} \]

计算过程是：

\[C_{1,1}=\max\{1+3,1+1\}=4\\ C_{1,2}=\max\{1+0,1+1\}=2\\ C_{2,1}=\max\{1+3,2+1\}=4\\ C_{1,1}=\max\{1+0,2+1\}=3 \]

同普通矩阵乘法一样，当矩阵 \(A\) 的列数与矩阵 \(B\) 的行数相同时，乘法才有意义。

（2）结合律

要将广义矩阵乘法应用到DP中，我们还要证明矩阵乘法满足结合律。

即

\[(A \times B )\times C=A \times (B \times C) \]

记 \(P=A\times B\)，有 \(P_{i,j}=\displaystyle\max_{k=1}^{n}({A_{i,k}+B_{k,j}})\) ；同理记 \(P'=B\times C\) ，有\(=P'_{i,j}=\displaystyle\max_{k=1}^{n}({B_{i,k}+C_{k,j}})\) 。

根据定义可以得到

\((P\times C)_{a,b}=\displaystyle\max_{p=1}^{n}({P_{a,p}+C_{p,b}})=\displaystyle\max_{p=1}^{n}({\displaystyle\max_{k=1}^{n}({A_{a,k}+B_{k,p}})+B_{p,b}})=\displaystyle\max_{p=1,k=1}^{n}({{A_{a,k}+B_{k,p}}+C_{p,b}})\)

\((A\times P')_{a,b}=\displaystyle\max_{p=1}^{n}({A_{a,p}+P'_{p,b}})=\displaystyle\max_{p=1}^{n}({\displaystyle\max_{k=1}^{n}({B_{p,k}+C_{k,b}})+A_{a,p}})=\displaystyle\max_{p=1,k=1}^{n}({{A_{a,p}+B_{p,k}}+C_{k,b}})\)

两式完全等价，所以满足结合律。不必理解证明过程，只要记住结论即可。

（3）交换律

不满足。所以一定要记住运算顺序。

（4）单位矩阵

如同数字的乘法有单位 \(1\)，广义矩阵乘法也有单位元，记为 \(I\) 。

\[I_{i,j}=\begin{cases} 0 &i=j\\ -\infty &otherwise \end{cases} \]

也就是说，单位矩阵的斜对角线为 \(0\)，其余为负无穷。很容易证明任何矩阵与之相乘都能得到其本身。

三、广义矩阵乘法的DP运用

回到最开始的问题，我们该如何在修改的情况下求出最长子段？

我们不妨将一开始的DP式写成矩阵的形式。同时，我们引入一个 \(g_i\) ，表示 \(\displaystyle\max_{j=1}^{i}f_j\) ，就是序列前 \(i\) 个值的答案。

于是有

\[f_i=\max\{f_{i-1}+a_{i},a_i\}\\ g_i=\max\{g_{i-1},f_i\}=\max\{g_{i-1},f_{i-1}+a_i,a_i\} \]

根据这个式子，我们可以构造矩阵 \(F_i=\begin{bmatrix}f_i\\g_i\\0 \end{bmatrix}\) 表示一组DP值。可以这样构造转移：

\[F_i=\begin{bmatrix} f_{i}\\ g_{i}\\ 0 \end{bmatrix}= \begin{bmatrix} a_i &-\infty& a_i\\ a_i&0&a_i\\ -\infty&-\infty&0 \end{bmatrix}\times \begin{bmatrix} f_{i-1}\\ g_{i-1}\\ 0 \end{bmatrix} \]

其实就是把转移式换了种写法对不对？很容易理解的。

发现矩阵式里有个矩阵 \(\begin{bmatrix}a_i &-\infty& a_i\\a_i&0&a_i\\-\infty&-\infty&0\end{bmatrix}\) ，仅和 \(i\) 有关，我们把它记为 \(A_i\) 。

可以有

\(F_n=A_nF_{n-1}=A_nA_{n-1}F_{n-2}=A_nA_{n-1}\cdots A_1F_0\)

初始状态 \(F_0=\begin{bmatrix} -\infty\\-\infty\\0\end{bmatrix}\)。（请注意一下矩阵相乘顺序）

聪明的同学发现，这不就是区间乘嘛！那么修改 \(a_i\) 就可以通过线段树来完成。我们用线段树维护区间信息的同时，也可以做到询问 \([l,r]\) 的最大子段。

这样做是对的吗？注意，我们刚刚已经证明了广义矩阵乘法满足结合律，所以完全可以维护区间矩阵乘。

这便是动态DP的核心思想。

这里提供这题的动态DP Code：

#include <bits/stdc++.h>
#define ll long long
using namespace std;

const ll inf=0x3f3f3f3f3f3f3f3fLL;
const int maxn=2e5+5;
ll n,a[maxn],m;

struct matrix{
	ll a[3][3];
	matrix(){
		memset(a,0xc0,sizeof a);
		for(int i=0;i<=2;i++){
			a[i][i]=0;
		}
	}
	void init(){
		memset(a,0,sizeof a);
	}
	void clear(){
		memset(a,0xc0,sizeof a);
	}
	matrix operator * (const matrix &x)const{
		matrix c;
		c.clear();
		for(int i=0;i<=2;i++){
			for(int j=0;j<=2;j++){
				for(int k=0;k<=2;k++){
					c.a[i][j]=max(c.a[i][j],a[i][k]+x.a[k][j]);
				}
			}
		}
		return c;
	}
};

struct seg_tree{
	int l,r;
	matrix sum;
}t[maxn*4];

void build(int o,int l,int r){
	t[o].l=l,t[o].r=r;
	if(l==r){
		matrix c;
		c.a[0][0]=c.a[0][2]=c.a[1][0]=c.a[1][2]=a[l];
		c.a[0][1]=c.a[2][0]=c.a[2][1]=-inf,c.a[1][1]=c.a[2][2]=0;
		t[o].sum=c;
		return;
	}
	int mid=(l+r)>>1;
	build(o<<1,l,mid);
	build(o<<1|1,mid+1,r);
	t[o].sum=t[o<<1|1].sum*t[o<<1].sum;
}

matrix query(int o,int l,int r){
	if(l<=t[o].l&&t[o].r<=r){
		return t[o].sum;
	}
	matrix c;
	if(r>=t[o<<1|1].l) c=c*query(o<<1|1,l,r);
	if(l<=t[o<<1].r) c=c*query(o<<1,l,r);
	return c;
}

void update(int o,int p,int k){
	if(t[o].l==t[o].r){
		t[o].sum.a[0][0]=t[o].sum.a[0][2]=t[o].sum.a[1][0]=t[o].sum.a[1][2]=k;
		return;
	}
	if(p<=t[o<<1].r) update(o<<1,p,k);
	else update(o<<1|1,p,k);
	t[o].sum=t[o<<1|1].sum*t[o<<1].sum;
}

int main(){
	scanf("%lld",&n);
	for(int i=1;i<=n;++i){
		scanf("%lld",&a[i]);
	}
	build(1,1,n);
	scanf("%lld",&m);
	while(m--){
		int op;
		scanf("%d",&op);int l,r;
		scanf("%d%d",&l,&r);
		if(op==0){
			update(1,l,r);
		}else{
			matrix ans;ans.init();
			ans.a[1][0]=ans.a[0][0]=-inf;
			ans=query(1,l,r)*ans;
			printf("%lld\n",ans.a[1][0]);
		}
	}
} 

四、树上问题

我们都知道树上最大独立集问题。如果不知道的话，就是“没有上司的舞会”带点权版。

现在我们有若干次修改点权，求每次修改后的最大独立集。

先不考虑修改，我们记 \(f_{u,0}\) 为不选择 \(u\) 的子树最大独立集；类似地，\(f_{u,1}\) 表示选择。

有

\[\begin{cases} f_{u,1}=a_i+\sum f_{v,0}\\ f_{u,0}=\sum \max(f_{v,0},f_{v,1}) \end{cases} \]

我们想要仿照上面的方法进行构造矩阵，但发现这个DP式有点臃肿。

为了满足修改要求，我们更改一下DP的转移，

我们定义 \(g_{u,1}\) 表示选择 \(u\) 这个点，只考虑 \(u\) 的轻子树的最大独立集；

\(g_{u,0}\) 表示不选择 \(u\) ，只考虑 \(u\) 的轻子树的最大独立集。

对 \(f\) 的定义不改变，假如我们已经求出了 \(g\)，那么转移式改为

\[\begin{cases} f_{u,1}=g_{u,1}+f_{v,0} \\f_{u,0}=g_{u,0}+\max(f_{v,0},f_{v,1}) \end{cases} \]

其中 \(v\) 是 \(u\) 的重子节点。

那么我们可以愉快地构造矩阵了：

\[F_u=\begin{bmatrix} f_{u,1}\\ f_{u,0} \end{bmatrix}= \begin{bmatrix} -\infty&g_{u,1}\\ g_{u,0}&g_{u,0} \end{bmatrix} \begin{bmatrix} f_{v,1}\\ f_{v,0} \end{bmatrix} \]

可以看到DP始终在重链上进行。我们惊喜地发现，系数矩阵 \(\begin{bmatrix}-\infty&g_{u,1}\\g_{u,0}&g_{u,0}\end{bmatrix}\) 只与 \(g\) 有关。可以这样考虑，对于每次修改，我们只需往上在重链里修改 \(g\) 值以及系数矩阵即可。问题来了，怎么修改 \(g\) 呢？

我们发现啊，由于 \(g\) 是只考虑轻子树的最大独立集，所以对于一个在重链上节点 \(u\) 的点权更新，只会影响到 \(g_{u,1}\) 以及往上跳时 \(fa_{top_u}\) 的信息。

我们还没说过怎么求 \(g\)，\(g\) 的大小只与轻儿子有关，

\[\begin{cases} g_{u,1}=a_u+\sum f_{v,0} \\g_{u,0}=\sum \max(f_{v,0},f_{v,1}) \end{cases} \]

其中 \(v\) 是轻儿子。当我们更新 \(u\) 时，我们能更新出从叶子节点到 \(top_u\) 的矩阵。利用这些矩阵，我们可以计算出 \(F_{top_u}=\begin{bmatrix} f_{top_u,1}\\f_{top_u,0}\end{bmatrix}\) ，并得出\(F_{top_u}\)的更新量。这个更新量就可以用于更新 \(g_{fa_{top_u}}\) 。

注意到

\(F_{u}=A_uF_v=A_uA_v\cdots F_{leaf}\)

而 \(F_{leaf}=A_{leaf}\begin{bmatrix}0\\0 \end{bmatrix}\)。

因为

\[\begin{bmatrix} -\infty &g_{u,1}\\ g_{u,0} &g_{u,0} \end{bmatrix} \times \begin{bmatrix} 0\\ 0 \end{bmatrix}= \begin{bmatrix} g_{u,1}\\ g_{u,0} \end{bmatrix} \]

当 \(u\) 为叶子，\(\begin{bmatrix}g_{u,1}\\g_{u,0}\end{bmatrix}=\begin{bmatrix}f_{u,1}\\f_{u,0}\end{bmatrix}\)

这个矩阵是 \(2\times 1\) 的，实际上可以看成是 \(2\times 2\) 矩阵的右半部分。

所以 \(F_u=A_uA_v\cdots A_{leaf}\begin{bmatrix}0\\0 \end{bmatrix}\)，实际上就是\(A_uA_v\cdots A_{leaf}\) 的右半部分。于是拿区间乘结果去更新 \(fa_{top_u}\) 即可。求最大独立集也就是求 \(F_1\) 。通过重链剖分与线段树实现，时间复杂度为 \(O(m\log^2n)\)。

这就是树上动态DP了。

这里提供P4719 【模板】"动态 DP"&动态树分治的Code:

#include <bits/stdc++.h>
#define ll long long
#define inf 0x3f3f3f3f3f3f3f3fLL
using namespace std;

const int maxn=1e5+5;
int fa[maxn],sz[maxn],son[maxn],top[maxn],tail[maxn],dfn[maxn],rnk[maxn],clk;
int a[maxn],n,m;
ll f[maxn][2];

struct node{
	int v,nxt;
}edge[maxn*2];int head[maxn],cntE;

struct matrix{
	ll a[2][2];
	matrix(){
		memset(a,0xc0,sizeof a);
		for(int i=0;i<=1;i++) a[i][i]=0;
	}
	void clear(bool k){
		if(k)memset(a,0xc0,sizeof a);
		else memset(a,0,sizeof a);
	}
	matrix operator * (const matrix &x) const{
		matrix c;c.clear(0);
		for(int i=0;i<2;i++)
			for(int j=0;j<2;j++)
				for(int k=0;k<2;k++)
					c.a[i][j]=max(c.a[i][j],a[i][k]+x.a[k][j]);
		return c;
	}
}g[maxn];

struct seg_tree{
	matrix sum;
	int l,r;
}t[maxn*4];

void init(int u,int p){
	fa[u]=p;sz[u]=1;
	f[u][0]=0;f[u][1]=a[u];
	for(int i=head[u];i;i=edge[i].nxt){
		int v=edge[i].v;
		if(v==p) continue;
		init(v,u);sz[u]+=sz[v];
		f[u][0]+=max(f[v][1],f[v][0]);
		f[u][1]+=f[v][0];
		if(!son[u]||sz[v]>sz[son[u]]) son[u]=v;
	}
}

void dfs(int u,int t){
	top[u]=t;
	dfn[u]=++clk;rnk[clk]=u;tail[t]=clk;
	g[u].a[0][1]=a[u];g[u].a[0][0]=-inf;
	if(!son[u]) return;
	dfs(son[u],t);
	for(int i=head[u];i;i=edge[i].nxt){
		int v=edge[i].v;
		if(v==fa[u]||v==son[u]) continue;
		g[u].a[0][1]+=f[v][0];
		g[u].a[1][1]+=max(f[v][0],f[v][1]);
		dfs(v,v);
	}
	g[u].a[1][0]=g[u].a[1][1];
}

void build(int o,int l,int r){
	t[o].l=l,t[o].r=r;
	if(l==r){
		t[o].sum=g[rnk[l]];
		return;
	}
	int mid=(l+r)>>1;build(o<<1,l,mid);build(o<<1|1,mid+1,r);
	t[o].sum=t[o<<1].sum*t[o<<1|1].sum;
}

matrix query(int o,int l,int r){
	if(l<=t[o].l&&t[o].r<=r){
		return t[o].sum;
	}
	matrix c;
	if(l<=t[o<<1].r) c=c*query(o<<1,l,r);
	if(r>=t[o<<1|1].l) c=c*query(o<<1|1,l,r);
	return c;
}

void modify(int o,int p){
	if(t[o].l==t[o].r){
		t[o].sum=g[rnk[p]];return;
	}
	if(p<=t[o<<1].r) modify(o<<1,p);
	if(p>=t[o<<1|1].l) modify(o<<1|1,p);
	t[o].sum=t[o<<1].sum*t[o<<1|1].sum;
}

void update(int u,int val){
	g[u].a[0][1]+=val-a[u];
	a[u]=val;
	while(u){
		matrix x=query(1,dfn[top[u]],tail[top[u]]);
		modify(1,dfn[u]);
		matrix y=query(1,dfn[top[u]],tail[top[u]]);
		u=fa[top[u]];
		g[u].a[1][1]+=max(y.a[0][1],y.a[1][1])-max(x.a[0][1],x.a[1][1]);
		g[u].a[0][1]+=y.a[1][1]-x.a[1][1];
		g[u].a[1][0]=g[u].a[1][1];
	}
}

void add(int u,int v){
	edge[++cntE]=(node){v,head[u]};
	head[u]=cntE;
}

int main(){
	scanf("%d%d",&n,&m);
	for(int i=1;i<=n;i++) scanf("%d",&a[i]);
	for(int i=1;i<=n-1;i++){
		int u,v;scanf("%d%d",&u,&v);
		add(u,v);add(v,u);
	}
	init(1,0);
	dfs(1,1);
	build(1,1,n);
	while(m--){
		int x,y;scanf("%d%d",&x,&y);
		update(x,y);
		matrix ans=query(1,1,tail[1]);
		printf("%lld\n",max(ans.a[1][1],ans.a[0][1]));
	}
}