【学习笔记】浅谈斜率优化 dp

· · 算法·理论

写一篇学习笔记,以防自己忘记。

引入

我们以 HDU 上的一道题为例。

Vjudge 提交入口

:::info[简易版题意] 现在有一个长度为 N 的序列 A,现在要将其分成若干个连续的组,定义一组的代价是这组内所有的 A 的和的平方再加上一个固定常数 M。求最小的总代价。

数据范围满足 0\le N\le500000,0\le M\le1000,要多测。 :::

暴力做法

显然的 dp,我们考虑这样设计状态:f_i 表示 ii+1 断开的情况下,1\sim i 的最小总费用。先定义数组 s 为数组 A 的前缀和。转移方程很好想,这里不阐述原因:

f_i=\min_{j=0}^{i-1}{f_j+(s_i-s_j)^2+M}

答案就是 f_N,时间复杂度显然是 O(n^2) 级别的。

斜率优化

这个式子不是很好的优化啊!咋办?

先竭尽所能的化简:

f_i=\min_{j=0}^{i-1}{f_j+s_i^2+s_j^2-2s_is_j}+M

注意到 s_i^2M 都是常数项,我们考虑将其转化为一次函数状物。

f_j+s_j^2=2s_is_j+f_i-s_i^2-M

我们印象中的一次函数长这样:

y=kx+b

说明一下其中的对应关系:

y=f_j+s_j^2,k=2s_i,x=s_j,b=f_i-s_i^2-M

注意到这里笔者特意将含 j 的都放到了 xy 当中,原因待会解释。

注意到这里 \min 消失了,我解释一下啥情况。就是我上面的将 x=s_j,y=f_j+s_j^2,可以将其转化为点。如果我们上面的一次函数穿过了这个点,由于斜率为 2s_i,我们可以直接解出 b,由于 -s_i^2-M 是常数项,也就可以算出 f_i 的大小。那这里我们需要 f_i 最小,也就是要让 b 最小。而 b 最小就是要让截距最小,这里放一张图。

注意到我上面话的意思就是要让那条直线向上平移最短的距离,碰到的第一个点就是能够让 b 取最小的点,这个点往往被我们叫做最优决策点

:::info[下凸壳的定义]

  1. 几何直观定义: 包围所有点的一条最底部的折线,使得没有任何一个点在这条折线的下方。
  2. 代数本质定义(写代码的核心): 从左到右,相邻两点连线的斜率严格单调递增(即 k_1 < k_2 < k_3 \dots)。

    一句话总结: 下凸壳就是斜率越来越大的底部折线。 :::

:::info[为什么最优的点一定出现在下凸壳上面] 假设有三个点 A, B, C,横坐标满足 x_A < x_B < x_C

如果中间的 B 点破坏了下凸壳,几何上意味着 B 点“凹进去了”,也就是它位于 AC 连线的上方。代数上则意味着,线段 AB 的斜率大于等于线段 BC 的斜率(即 k_{AB} \ge k_{BC})。

我们在求最优解时,相当于用一条斜率为 k 的直线从下往上平移,寻找第一个碰到的点(即最小截距):

  1. 如果直线的斜率 k 较小(即 k < k_{AB}),直线在平移碰到 B 之前,一定会先碰到左边的 A
  2. 如果直线的斜率 k 较大(即 k \ge k_{BC}),直线在平移碰到 B 之前,一定会先碰到右边的 C

因为 k_{AB} \ge k_{BC},所以无论直线的斜率 k 取什么值,它在碰到 B 之前,必然已经先碰到了 AC

因此,位于连线上方的 B 点在任何情况下都不可能成为最优解,是一个废点。把所有这样的废点剔除后,剩下的点自然就构成了斜率严格单调递增的下凸壳。 :::

我们现在要维护这个下凸壳。注意到我们的 s 由于是前缀和,因此随 i 的增大,s_i 也在增大,所以这条线的斜率是不断增大的。我们再来看一下那个点有什么特征,注意到该点左边的斜率全部小于 2s_i,右边的斜率全部大于 2s_i,所以维护的方法显而易见了:我们使用一个单调队列,由于 k 是不断增加的,所以这一条曲线的左下角如果斜率小于 k 就直接删除,最后留下的队头就是我们心心念念的 j
找到 j 以后计算 f_i 的方式有两种。一种是使用我们先前的暴力转移式来转移,一种是将我们刚才的一次函数带入算出 b 并算出 f_i,这两种方法都是可行的。

这里其实还没结束。我们算出 f_i,还需要往这个图中加入 (s_i,f_i+s_i^2),由于这个横纵坐标显然都是单调递增的,因此这个点会出现在右上角,并要将其与上一个点连边。

我们现在完全盯着上面那个图:新加的 i-1 点要与 t 连边,但是,如果 t-1\rightarrow t 的斜率大于等于 t\rightarrow i-1 的斜率,这就意味着 t 点凹进去了,下凸壳结构被破坏了!所以我们需要用一个 while 循环判断,如果破坏了结构就将队尾弹出,直到满足单调递增为止,然后再将 i 塞进队列。

这一块确实有点抽象,建议没看懂的同学自己在脑海里面多过几遍。

代码特别好写了: ::::success[code] :::warning[关于斜率的计算部分]{open} 如果我们现在要计算 AB 的斜率,但是 x_A=x_B 要怎么办?我们都知道对于满足 x_A\le x_BAB 的斜率的计算方式为 \frac{\Delta y}{\Delta x},也就是 \frac{y_B-y_A}{x_B-x_A},但是当 x_A=x_B 时,分母为 0。比较简单的处理方式有两种,一种是使用交叉相乘大法,原本的判断条件:\frac{y_B - y_A}{x_B - x_A} \ge \frac{y_C - y_B}{x_C - x_B} 改写为:(y_B - y_A) \times (x_C - x_B) \ge (y_C - y_B) \times (x_B - x_A)。另一种是根据 y_Ay_B 的关系给它赋值个极值,意味着当 x_A=x_B 时的斜率是最大/最小的。

所以我们在写题时如果意识到这道题目中的横坐标可能相同,那么记得对其作出处理!不然可能会喜提 RE。 ::: :::warning[一个小细节]{open} 注意到我的首先单调队列初始情况下 h=t=0,原本应该是 h=1,t=0 的,为什么这样做呢?其实这项等于初始情况下塞进去了一个 0,注意到我们暴力的转移式 j 可以等于 0,因此要提前往里面塞一个 0,方便转移。后面的代码都有这个细节。 :::

#include <bits/stdc++.h>
using namespace std;
#define int long long
const int N=5e5+5;
int n,m,f[N],q[N],s[N],head,tail,x;
double slope(int x,int y,int xx,int yy){return x==xx?(yy>=y?1e18:-1e18):1.0*(yy-y)/(xx-x);}
int get_x(int x){return s[x];}
int get_y(int x){return f[x]+s[x]*s[x];}
signed main(){
    while(cin>>n>>m){
        head=tail=0;
        for(int i = 1;i<=n;i++)
            cin>>x,s[i]=s[i-1]+x;
        for(int i = 1;i<=n;i++){
            int k=2*s[i];
            while(head<tail&&slope(get_x(q[head]),get_y(q[head]),get_x(q[head+1]),get_y(q[head+1]))<k) head++;
            f[i]=f[q[head]]+(s[i]-s[q[head]])*(s[i]-s[q[head]])+m;
            while(head<tail&&slope(get_x(q[tail-1]),get_y(q[tail-1]),get_x(q[tail]),get_y(q[tail]))>=slope(get_x(q[tail]),get_y(q[tail]),get_x(i),get_y(i))) tail--;
            q[++tail]=i;
        }
        cout<<f[n]<<'\n';
    }
    return 0;
}

::::

时间复杂度显然是 O(n)

然后我们发现斜率优化的大部分题就是列出式子,然后再转化为一次函数的形式,剩下的就是模版了。

拓展

P5785 [SDOI2012] 任务安排

我们先来思考一下这道题目怎么做。

首先先列出来暴力的转移式长什么样子(定义 c 数组为费用的前缀和,定义 t 数组为时间的前缀和):

f_i=\min_{j=0}^{i-1}{f_j+s(c_n-c_j)+t_i(c_i-c_j)}

这个式子就是先将这一次分组的 s 对后面产生的花费提前计算,并加上这一组本身需要的花费。

化简得:

f_i=f_j+s \cdot c_n-s \cdot c_j+t_i \cdot c_i-t_i \cdot c_j

移项得:

f_j=s \cdot c_j+t_i \cdot c_j+f_i-s \cdot c_n-t_i \cdot c_i\\ f_j=(s+t_i) \cdot c_j+f_i-s \cdot c_n-t_i \cdot c_i

其中 y=f_j,k=s+t_i,x=c_j,b=f_i-s \cdot c_n-t_i \cdot c_i

接下来就是模版的事了。 ::::success[code] :::warning[该代码疑似有被 hack 的风险]{open} 注意到该代码中计算斜率的函数并没有判断 x 相等,但是在题目的意思下,花费的费用可能为 0,也就是前缀和数组 c 中可能有相同的数值。但是这份代码通过了此题,并且由于笔者太懒,并没有对此作出修改。 :::

#include <bits/stdc++.h>
using namespace std;
#define int long long
const int N=3e5+5;
int n,s,t[N],c[N],x[N],y[N],f[N],q[N],head,tail;
double slope(int x,int y,int xx,int yy){return 1.0*(yy-y)/(xx-x);}
signed main(){
    cin>>n>>s;
    for(int i = 1;i<=n;i++)
        cin>>x[i]>>y[i],t[i]=t[i-1]+x[i],c[i]=c[i-1]+y[i];
    for(int i = 1;i<=n;i++){
        int k=s+t[i];
        while(head<tail&&slope(c[q[head]],f[q[head]],c[q[head+1]],f[q[head+1]])<=k) head++;
        f[i]=f[q[head]]+s*(c[n]-c[q[head]])+t[i]*(c[i]-c[q[head]]);
        while(head<tail&&slope(c[q[tail-1]],f[q[tail-1]],c[q[tail]],f[q[tail]])>=slope(c[q[tail]],f[q[tail]],c[i],f[i])) tail--;
        q[++tail]=i;
    }
    cout<<f[n];
    return 0;
}

:::: 好的现在我们再来看一下当时间不一定为正数咋做 (吐槽一下,时间为负数这合理吗)。时间不一定为正数这意味着我们的前缀数组 t 不一定是单调递增的,因此被破坏单调性的有我们的斜率 k=s+t_i,还有纵坐标 y=f_j。我们先来看一下 y=f_j 的单调性被破坏有影响吗,这就意味着可能出现斜率为负数的线,但是无论它斜率是不是负数,只要满足 t-1\rightarrow t 的斜率大于 t\rightarrow i-1 就会破坏下凸壳结构,都要弹出。所以 y 的单调性被破坏是没影响的。但是 k 的单调性被破坏,就意味着我们不能再使用单调队列 O(n) 来维护了,因为如果提前弹出会导致后面斜率更小的一次函数找不到最优决策点,因此我们不能直接使用单调队列来维护。注意到我们寻找到的最优决策点左部分的斜率均满足 <k,右部分的斜率均满足 >k,因此我们可以考虑直接使用二分答案来维护。然后其实这个单调队列的队头我们是永远不会去管它的,所以这个单调队列或许在这里其实是一个单调栈。

时间复杂度显然的,每一次都要一个二分,O(n\log n)

:::success[code]

#include <bits/stdc++.h>
using namespace std;
#define int long long
const int N=3e5+5;
int n,s,x,y,c[N],t[N],q[N],f[N],head,tail;
double slope(int x,int y,int xx,int yy){return x==xx?(yy>=y?1e18:-1e18):1.0*(yy-y)/(xx-x);}
signed main(){
    cin>>n>>s;
    for(int i = 1;i<=n;i++)
        cin>>x>>y,t[i]=t[i-1]+x,c[i]=c[i-1]+y;
    for(int i = 1;i<=n;i++){
        int l=head,r=tail-1,res=tail,k=s+t[i];
        while(l<=r){
            int mid=l+r>>1;
            if(slope(c[q[mid]],f[q[mid]],c[q[mid+1]],f[q[mid+1]])>k) r=mid-1,res=mid;
            else l=mid+1;
        }
        f[i]=f[q[res]]+s*(c[n]-c[q[res]])+t[i]*(c[i]-c[q[res]]);
        while(head<tail&&slope(c[q[tail-1]],f[q[tail-1]],c[q[tail]],f[q[tail]])>=slope(c[q[tail]],f[q[tail]],c[i],f[i])) tail--;
        q[++tail]=i;
    }
    cout<<f[n];
    return 0;
}

:::

其实还有一种被破坏单调性的可能,就是横坐标 x 的单调性破坏。这意味着我们不能单单只维护队尾,这时候我们就不能再用单调队列这种数据结构了,要换更高级的平衡树/cdq 分治/李超线段树来维护了,这一部分由于笔者还没学,等学了以后会补的!qwq

撰写不易,能给个赞再走吗 qwq