斜率优化小结

/ 0评 /

在一维线性动态规划里,我们经常见到形如 f(i)=min/max(f(j)+c_i)(其实可以看成 f(i)=min/max(a_i+b_j) )这样形式的转移方程。朴素做法是 \Theta (N^2)。显然这样的形式可以用单调队列维护,优化到 \Theta (N)。但是如果遇到 f(i)=min/max(a_i\ast b_j+c_i+d_j) 这样的,似乎用单调队列就难以维护了……因为 a_i\ast b_j 既与 i 有关,又与 j 有关。这时候就要用到另一种优化方式:斜率优化。

本文含有大量 KaTeX 公式,请确保浏览器资磁。
可能是我所有博客里公式最多的一篇,在编辑本文的时候,编辑器已经卡得不成样子了……

概述

下面我们以 f(i)=a_i\ast b_j+c_i+d_j 这个状态转移方程为例讲讲通用做法(如果不想看可以直接看例题)。

f(i)=a_i\ast b_j+c_i+d_j 这个状态转移方程,既和 i 有关又和 j 有关。简单来说,这个方程通过适当的转换可以看成一条直线,一般方法是令 x=b_j,y=d_j,通过变形就可以得到:

y=-a_i x+(f(i)-c_i)

可以看到这就是 y=kx+b 的形式,也就是说这个等式表示的是一条斜率为 -a_i 并且截距为 f(i)-c 的直线。既然要 f(i) 尽量小,自然要使这条直线截距尽量小。现在斜率确定了,那么我们要找之前的可以使截距最小的 (b_j,d_j)。通过一个单调队列维护下凸壳,我们可以在 \Theta (1) 时间内做到这一点。具体方法下面会讲。

例题

HNOI2008 玩具装箱TOY

洛谷题目链接:HNOI2008 玩具装箱TOY

题目描述

P教授要去看奥运,但是他舍不下他的玩具,于是他决定把所有的玩具运到北京。他使用自己的压缩器进行压缩,其可以将任意物品变成一堆,再放到一种特殊的一维容器中。P教授有编号为 1\cdots NN 件玩具,第 i 件玩具经过压缩后变成一维长度为 C_i .为了方便整理,P教授要求在一个一维容器中的玩具编号是连续的。同时如果一个一维容器中有多个玩具,那么两件玩具之间要加入一个单位长度的填充物,形式地说如果将第 i 件玩具到第 j 个玩具放到一个容器中,那么容器的长度将为 x=j-i+\sum\limits_{k=i}^{j}C_k 制作容器的费用与容器的长度有关,根据教授研究,如果容器长度为 x ,其制作费用为 (X-L)^2 .其中 L 是一个常量。P教授不关心容器的数目,他可以制作出任意长度的容器,甚至超过 L 。但他希望费用最小.。

数据范围:1\leqslant N\leqslant 50000,1\leqslant L,C_i\leqslant 10^7

分析

这是一道斜率优化的入门题,非常典型。

首先根据题目可以推出一维动态规划的状态转移:

f(i)=f(j)+(i-j+\sum_{k=j}^{i}(C_k)-L-1)^2

到这里可以 \Theta(N^2) 做。接下来前缀和累计一下,得到:

f(i)=f(j)+(i+sum_i-j-sum_{j-1}-L-1)^2

方便起见,令 A_i=i+sum_i,B_i=i+sum_{i-1}+L+1,原式又变成:

f(i)=f(j)+(A_i-B_j)^2=f(j)+2A_iB_j+A_i^2+B_j^2 \\
2A_iB_j+A_i^2+f(i)=f(j)+B_j^2

现在我们令 x=B_j,y=f(j)+B_j^2,那么得到:

y=2A_ix+(f(i)+A_i^2)

可以看到这就是 y=kx+b 的直线表达式形式。我们可以把它看做一条经过点 (B_j,f(j)+B_j^2) 并且斜率为 2A_i 的直线。那么现在我们要让 f(i) 尽量小,也就是这条直线在 y 轴上的截距尽量小。
想象在平面直角坐标系里这样斜率固定的这样一条直线,我们只需要一直平移平移平移,第一个与这条直线相交的 (B_j,f(j)+B_j^2) 的点就是满足截距最小的点。

在上图中,A 点既是当前最合适的点。但是怎么才能 \Theta (1) 找到这个点呢?
其实可以用单调队列维护一个下凸壳,如下图:

我们只对这个凸壳上的点感兴趣。那么具体怎么维护这样一个凸壳呢?我们维护一个单调队列,对于当前的 i,我们先把 head 开始不符合的点都删掉,删到 head 为最优为止,然后通过 head 算出当前答案;然后把不符合条件的 tail 都搞掉,再加入当前的点。很类似计算几何里求凸包。
具体实现,假设 slope(i,j) 代表 i,j 点连边的斜率:

删除头部的过程中,如果 slope(que_{head},que_{head+1})<2\ast A_i 则舍弃头,head++;
如果不能推头了,用当前 head 修正 f(i)
删尾部,如果 slope(que_{tail-1},que_{tail}< slope(que_{tail-1},i)) 就说明尾部没有当前这个点优秀,删尾部,删到不能删为止。

OK 了。复杂度从 \Theta (N^2) 降到了 \Theta (N)

代码

#include<cstdio>
#include<cstring>
#include<iostream>
#include<vector>
#define CLEAR(_) memset(_,0,sizeof(_))
#define sqr(_) ((_)*(_))
using namespace std;

const int maxn=50005;
int n;
long long f[maxn],sum[maxn],a[maxn],b[maxn],L;

struct Point{
    double x,y;
    int id;
};

int head=0,tail=0;
Point que[maxn];

inline int read(){
    int ret=0,f=1;char ch=getchar();
    while (ch<'0'||ch>'9'){if (ch=='-') f=-1;ch=getchar();}
    while (ch>='0'&&ch<='9') ret=ret*10+ch-'0',ch=getchar();
    return ret*f;
}

inline double GetK(Point p1,Point p2){
    return (double)(p1.y-p2.y)/(p1.x-p2.x);
}

int main(){
    n=read();L=read();
    b[0]=L+1;
    for (int i=1;i<=n;i++){
        int x=read();
        sum[i]=sum[i-1]+(long long)x;
        a[i]=sum[i]+i;
        b[i]=sum[i]+i+(long long)L+1.0;
    }
    head=tail=1;que[1]=(Point){b[0],sqr(b[0]),0};
    for (int i=1;i<=n;i++){
        while ((head<tail)&&(GetK(que[head],que[head+1])<2.0*(double)a[i])) head++;
        int j=que[head].id;
        f[i]=f[j]+sqr(a[i]-b[j]);
        Point now=(Point){(double)b[i],(double)f[i]+sqr(b[i]),i};
        while ((head<tail)&&(GetK(now,que[tail-1])<GetK(que[tail-1],que[tail]))) tail--;
        que[++tail]=now;
    }
    printf("%lld\n",f[n]);
    return 0;
}

ZJOI2007 仓库建设

洛谷题目链接:ZJOI2007 仓库建设

L公司有N个工厂,由高到底分布在一座山上。

工厂1在山顶,工厂N在山脚。 由于这座山处于高原内陆地区(干燥少雨),L公司一般把产品直接堆放在露天,以节省费用。

突然有一天,L公司的总裁L先生接到气象部门的电话,被告知三天之后将有一场暴雨,于是L先生决定紧急在某些工厂建立一些仓库以免产品被淋坏。

由于地形的不同,在不同工厂建立仓库的费用可能是不同的。第i个工厂目前已有成品Pi件,在第i个工厂位置建立仓库的费用是Ci。

对于没有建立仓库的工厂,其产品应被运往其他的仓库进行储藏,而由于L公司产品的对外销售处设置在山脚的工厂N,故产品只能往山下运(即只能运往编号更大的工厂的仓库),当然运送产品也是需要费用的,假设一件产品运送1个单位距离的费用是1。

假设建立的仓库容量都都是足够大的,可以容下所有的产品。你将得到以下数据:

请你帮助L公司寻找一个仓库建设的方案,使得总的费用(建造费用+运输费用)最小。

分析

显然是动态规划的思想,根据题意不难推出这个转移方程:

f(i)=f(j)+c_i+\sum_{k=j+1}^{i} P_k\ast dist(k,i)

一样的套路。
我们先定义几个量:

因为 dist(k,i)=d_i-d_{k-1},我们得到:

\sum_{k=j+1}^{i} P_k\ast dist(k,i)=\sum_{k=j+1}^{i} P_k(d_i-d_{k-1})=\sum_{k=j+1}^{i} P_kd_i-P_kd_{k-1}

展开得到:

d_i(P_{j+1}P_{j+2}\dots P_{i})-(P_{j+1}d_j+P_{j+2}d_{j+1}+\dots+P_i d_{i-1})

前缀和派上用场了,原式可以变形为:

f(i)=f(j)+c_i+d_i(sump_i-sump_j)-(summ_i-summ_j)

整理得到:

f(i)=c_i+d_i sump_i-summ_i-d_i sump_j+summ_j+f(j)

为了方便表述,我们令 A_i=c_i+d_i sump_i-summ_i,B_i=summ_i+f(i)

f(i)=A_i+B_j-d_i sump_j

然后令 x=sump_j,y=B_j 那么原式化为:

f(i)=A_i+y-d_ix \iff y=d_ix+f(i)-A_i

大功告成了!又是和上面那题一样的题目了。

代码

#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#include<vector>
using namespace std;
const int maxn=1000005;

int n;
long long d[maxn],c[maxn],p[maxn];
long long f[maxn],sump[maxn],summ[maxn];
long long a[maxn],b[maxn];
int head=0,tail=0,que[maxn];

inline int read(){
    int ret=0,f=1;char ch=getchar();
    while (ch<'0'||ch>'9'){if (ch=='-') f=-1;ch=getchar();}
    while (ch>='0'&&ch<='9') ret=ret*10+ch-'0',ch=getchar();
    return ret*f;
}

inline double X(int i){return (double)sump[i];}
inline double Y(int i){return (double)b[i];}

inline double slope(int i,int j){
    return (double)(Y(i)-Y(j))/(X(i)-X(j));
}

int main(){
    n=read();
    for (int i=1;i<=n;i++){
        d[i]=read();p[i]=read();c[i]=read();
        sump[i]=sump[i-1]+p[i];
        summ[i]=summ[i-1]+p[i]*d[i];
    }
    for (int i=1;i<=n;i++) a[i]=c[i]+d[i]*sump[i]-summ[i],f[i]=c[i]+d[i]*sump[i]-summ[i];
    head=tail=1;
    for (int i=1;i<=n;i++){
        while ((head<tail)&&(slope(que[head],que[head+1])<d[i])) head++;
        f[i]=min(f[i],a[i]+b[que[head]]-d[i]*sump[que[head]]);
        b[i]=summ[i]+f[i];
        while ((head<tail)&&(slope(que[tail-1],que[tail])>slope(que[tail-1],i))) tail--;
        que[++tail]=i;
    }
    printf("%lld\n",f[n]);
    return 0;
}

知识共享许可协议 本文章采用 知识共享署名-非商业性使用-相同方式共享 4.0 国际许可协议 进行许可。
欢迎转载,如有错误欢迎指出。
本文链接:https://skywt.cn/posts/slope/


发表评论

电子邮件地址不会被公开。 必填项已用*标注