斜率优化DP总结

时间:2024-01-27 14:46:04

前言

(本文中的图片都由\(WPS\)出品)
\(DP\)\(OI\) 中重要的一部分
一般来说,因为 \(DP\) 会把之前的结果保存下来,所以时间复杂度还是比较优秀的
但是在某些情况下,时间复杂度仍然超出了题目的限制,这是我们就要考虑对其进行优化
\(DP\) 的优化一般从状态、决策、转移三个方面去考虑
而斜率优化则是对决策进行优化的一种方法
它适用于类似 \(f[i]=min/max(a[i] \times b[j]+c[i]+d[j])\) 的方程

例题(洛谷 P4072 [SDOI2016]征途 )

题目描述

\(Pine\) 开始了从 \(S\) 地到 \(T\) 地的征途。

\(S\)地到\(T\)地的路可以划分成 \(n\) 段,相邻两段路的分界点设有休息站。

\(Pine\)计划用\(m\)天到达\(T\)地。除第\(m\)天外,每一天晚上\(Pine\)都必须在休息站过夜。所以,一段路必须在同一天中走完。

\(Pine\)希望每一天走的路长度尽可能相近,所以他希望每一天走的路的长度的方差尽可能小。

帮助\(Pine\)求出最小方差是多少。

设方差是\(v\),可以证明,\(v\times m^2\)是一个整数。为了避免精度误差,输出结果时输出\(v\times m^2\)

输入格式

第一行两个数 \(n\)\(m\)

第二行 \(n\) 个数,表示 \(n\) 段路的长度

输出格式

一个数,最小方差乘以 \(m^2\) 后的值

输入输出样例

输入 #1

5 2
1 2 5 8 6

输出 #1

36

说明/提示

对于 \(30\%\) 的数据,\(1 \le n \le 10\)

对于 \(60\%\) 的数据,\(1 \le n \le 100\)

对于 \(100\%\) 的数据,\(1 \le n \le 3000\)

保证从 \(S\)\(T\) 的总路程不超过 \(30000\)

分析

要对一个状态转移方程进行优化,首先要把最朴素的方程写出来

在本题中,稍加推导即可写出时间复杂度为 \(O(m \times n^2)\)的状态转移方程

\[s^2 \times m^2=\frac{(v_1-\overline{v})^2+(v_2-\overline{v})^2+...+(v_m-\overline{v})^2}m \times m^2 \]

\[s^2 \times m^2=((v_1-\overline{v})^2+(v_2-\overline{v})^2+...+(v_m-\overline{v})^2) \times m \]

\[s^2 \times m^2=m\times \sum_{i=1}^mv_i^2+m^2 \times \overline {v}^2-2\times m \times \overline {v} \times sum[n] \]

又因为 \(\overline{v}=\frac{sum[n]}{m}\)

所以

\[s^2\times m^2=m\times \sum_{i=1}^mv_i^2-sum[n]\times sum[n] \]

后面的值是固定的,所以我们只需要让前面的值最小化即可

我们设\(f[i][j]\)为前\(i\)天分成\(j\)段所得到的最小值

那么就有

\[f[i][k]=\min(f[j][k-1]+(sum[i]-sum[j])^2) \]

展开就有

\[f[i][k]=f[j][k-1]+sum[i]^2+sum[j]^2-2\times sum[i] \times sum[j] \]

接下来就是本文的重点:如何用斜率优化这类方程

首先,你需要掌握一次函数 \(y=kx+b\) 的图像和性质

这应该问题不大

下面我们就要对方程进行移项,使其变成易于优化的形式

\[f[j][k-1]+sum[j]^2=2\times sum[i] \times sum[j]+f[i][k]-sum[i]^2 \]

我们发现,这和一次函数的解析式完全吻合

我们把\(f[j][k-1]+sum[j]^2\)看成\(y\)

\(2 \times sum[i]\)看成\(k\)

\(sum[j]\)看成\(x\)

\(f[i][k]-sum[i]^2\)看成\(b\)

这样,对于每一个\(i\)来说,直线的\(k\)是确定的

我们要使\(f[i][k]\)最小,也就是要使\(b\)最小

我们可以把所有的\(j\)想象成空间坐标为\((sum[j],f[j][k-1]+sum[j]^2)\)中的点

知道了斜率,知道了直线上的点,那么这条直线就确定了

那么我们考虑什么样的点使直线的\(b\)最小

直线\(l\)是我们要移动的直线,平面中的点是可以转移的\(j\)

我们会发现当当前点和后一个点形成的直线的斜率恰好大于直线\(l\)的斜率时,由当前点转移决策是最优的

在这里要特别强调一下:本题中斜率\(k\) 和 横坐标 \(x\) 均为单调递增的,对于 \(k\)\(x\)不单调递增的情况,处理方式不同

这就是代码里面的

while(head<tail && xl(q[head],q[head+1])<2*sum[j]) head++;

我们再去考虑什么样的点肯定不会对结果产生贡献

上面的图中\(2\)号节点是无论如何也不会更新其它节点的

因为\(1\)号节点或\(3\)号节点总会比它更优

这就是代码里的

while(head<tail && xl(j,q[tail-1])>=xl(j,q[tail])) tail--;

整个过程就相当于维护了一个下凸包

在求斜率的函数中,我们要判掉 \(x\) 相等的情况,在某些时候,还要判掉 \(y\) 相等的情况

double X(int id){
	return (double)sum[id];
}
double Y(int id){
	return (double)(g[id]+sum[id]*sum[id]);
}
double xl(int i,int j){
	if(std::fabs(X(i)-X(j))<eps){
		if(std::fabs(Y(i)-Y(j))<eps) return 0;
		else if(Y(i)>Y(j)) return 1e18;
		else return -1e18;
	}
	return (Y(i)-Y(j))/(X(i)-X(j));
}

拓展一:斜率不单调但x单调

如果斜率不是单调递增,我们就不能从前面清空队列直接转移

比如上面这幅图如果在遇到直线 \(m\) 时一直从前清空队列的话那么就会把\(3\)号决策点弹出队列

但是如果之后遇到一个斜率比较小的直线\(l\)那么就不能转移到最优解

典型的例题是 洛谷P5785 [SDOI2012]任务安排

朴素的状态转移方程为

$ f[i] = f[j] + (sumc[i] - sumc[j]) * sumt[i] + s * (sumc[n] - sumc[j]);$

在这一道题中,作为斜率的 \(sumt\) 不再单调

但是 \(x\) 之仍然是单调的

所以我们可以用维护一个斜率单调的队列

每次在队列中二分答案

值得一提的是,出题人精心准备了卡精度的数据,所以我们要把二分时的除法改为乘法

代码

#include <cstdio>
#define rg register
inline int read() {
    rg int x = 0, fh = 1;
    rg char ch = getchar();
    while (ch < '0' || ch > '9') {
        if (ch == '-')
            fh = -1;
        ch = getchar();
    }
    while (ch >= '0' && ch <= '9') {
        x = (x << 1) + (x << 3) + (ch ^ 48);
        ch = getchar();
    }
    return x * fh;
}
typedef long long ll;
const int maxn = 1e6 + 5;
int t[maxn], c[maxn], n, s, q[maxn], head, tail;
ll f[maxn], sumt[maxn], sumc[maxn];
double Y(int i) { return (double)(f[i] - s * sumc[i]); }
double X(int i) { return (double)sumc[i]; }
double xl(int i, int j) {
    if (X(i) == X(j)) {
        if (Y(i) > Y(j))
            return 1e18;
        else
            return -1e18;
    }
    return (double)(Y(i) - Y(j)) / (X(i) - X(j));
}
int ef(double now) {
    int l = head, r = tail, mids;
    while (l < r) {
        mids = (l + r) >> 1;
        if ((X(q[mids]) > X(q[mids + 1]) &&
             Y(q[mids]) - Y(q[mids + 1]) < now * (X(q[mids]) - X(q[mids + 1]))) ||
            (X(q[mids]) < X(q[mids + 1]) &&
             (Y(q[mids + 1]) - Y(q[mids]) < now * (X(q[mids + 1]) - X(q[mids])))))
            l = mids + 1;
        else
            r = mids;
    }
    return q[l];
}
int main() {
    n = read(), s = read();
    for (rg int i = 1; i <= n; i++) {
        t[i] = read();
        c[i] = read();
        sumt[i] = sumt[i - 1] + t[i];
        sumc[i] = sumc[i - 1] + c[i];
    }
    head = tail = 1;
    for (rg int i = 1; i <= n; i++) {
        rg int wz = ef((double)(sumt[i]));
        f[i] = f[wz] + (sumc[i] - sumc[wz]) * sumt[i] + s * (sumc[n] - sumc[wz]);
        while (head < tail && xl(i, q[tail - 1]) >= xl(i, q[tail])) tail--;
        q[++tail] = i;
    }
    printf("%lld\n", f[n]);
    return 0;
}

扩展二、x不单调

没有找到 \(x\) 不单调但是 \(k\) 单调的题

但是却有一道 \(x\) 不单调 \(k\) 也不单调的题洛谷P4655 [CEOI2017]Building Bridges

对于这道题,我们同样可以写出最朴素的方程

\(f[i]f[j]+(h[i]-h[j])*(h[i]-h[j])+sum[i-1]-sum[j]\)

神奇的一点是 \(h\) 数组既作为直线的斜率又作为 \(x\)

而且 \(h\) 并不单调

这是我们就不能再用单调队列去维护,因为凸包的形状在不断改变

本题可以用李超线段树或者平衡树动态维护凸包解决

但是还有一种 \(CDQ\) 分治离线处理的方法

直接人为地排出单调性,像普通单调队列那样维护就可以了

如果左侧斜率递增,并且左侧编号小于右侧,那么可以通过单调队列维护左侧的凸包来更新右侧答案

并且这样一定能够遍历出每个节点的所有决策点

代码

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cstring>
#include<cmath>
#define rg register
inline int read(){
	rg int x=0,fh=1;
	rg char ch=getchar();
	while(ch<'0' || ch>'9'){
		if(ch=='-') fh=-1;
		ch=getchar();
	}
	while(ch>='0' && ch<='9'){
		x=(x<<1)+(x<<3)+(ch^48);
		ch=getchar();
	}
	return x*fh;
}
const int maxn=1e5+5;
const double eps=1e-6;
typedef long long ll;
int n,h[maxn],w[maxn];
ll f[maxn],sum[maxn];
double X(int id){
	return (double)h[id];
}
double Y(int id){
	return (double)(f[id]+1LL*h[id]*h[id]-sum[id]);
}
double xl(int i,int j){
	if(std::fabs(X(i)-X(j))<eps){
		if(std::fabs(Y(i)-Y(j))<eps) return 0;
		else if(Y(i)>Y(j)) return 1e18;
		else return -1e18;
	} else {
		return (Y(i)-Y(j))/(X(i)-X(j));
	}
}
bool cmp(int aa,int bb){
	return h[aa]<h[bb];
}
int tmp[maxn],p[maxn],q[maxn];
void solve(int l,int r){
	if(l==r) return;
	rg int mids=(l+r)>>1,head=l-1,tail=mids;
	for(rg int i=l;i<=r;i++){
		if(p[i]<=mids) tmp[++head]=p[i];
		else tmp[++tail]=p[i];
	}
	for(rg int i=l;i<=r;i++){
		p[i]=tmp[i];
	}
	solve(l,mids);
	head=1,tail=0;
	for(rg int i=l;i<=mids;i++){
		while(head<tail && xl(p[i],q[tail-1])>=xl(p[i],q[tail])) tail--;
		q[++tail]=p[i];
	}
	for(rg int i=mids+1;i<=r;i++){
		while(head<tail && xl(q[head+1],q[head])<=2.0*h[p[i]]) head++;
		f[p[i]]=std::min(f[p[i]],f[q[head]]+1LL*(h[p[i]]-h[q[head]])*(h[p[i]]-h[q[head]])+sum[p[i]-1]-sum[q[head]]);
	}
	solve(mids+1,r);
	head=l,tail=mids+1;
	for(rg int i=l;i<=r;i++){
		if(tail>r || (head<=mids && h[p[head]]<h[p[tail]])) tmp[i]=p[head++];
		else tmp[i]=p[tail++];
	}
	for(rg int i=l;i<=r;i++){
		p[i]=tmp[i];
	}
}
int main(){
	memset(f,0x3f,sizeof(f));
	n=read();
	for(rg int i=1;i<=n;i++){
		h[i]=read();
		p[i]=i;
	}
	for(rg int i=1;i<=n;i++){
		w[i]=read();
		sum[i]=sum[i-1]+w[i];
	}
	std::sort(p+1,p+1+n,cmp);
	f[1]=0;
	solve(1,n);
	printf("%lld\n",f[n]);
	return 0;
}