【一路走下去的斜率优化动态规划】

·随着网上众多OIer的步子,大米饼便静静地做了以下题目。

 

·首先列出大米饼的码风(代码风格):

①for循环被转化为Go循环和Ro循环分别表示升序和降序。②对于维护DP的单调队列,两个指针常用 Head和Tail两条。③对斜率优化一类题目的坐标点的宏定义X(i)Y(i),便于理解同时使用double Rate函数计算两点直线斜率。

 

[1]玩具装箱(详细阐述) 【LINK】

步骤一:

     列出DP方程式:设f[i]表示分组完前i件物品的最小花费,为方便计算,

设sum[i]表示是前i件物品的长度和。

     f[i]=Min(f[j]+(sum[i]-sum[j]+i-j-L-1)2)   [0<=j<i]

步骤二:

      对于范围内的j,进行式子变形,获得直线解析式:(先让L++)

      f[i]=f[j]+[(sum[i]+i)-(sum[j]+j)-L]2------->令s[k]=sum[k]+k

      f[i]=f[j]+(s[i]-s[j]-L)2----------------->整体法展开

      f[i]=f[j]+s[i]2+(s[j]+L)2-2*s[i]*(s[j]+L)---->移项

      f[i]+2*s[i]*(s[j]+L)=f[j]+s[i]2+(s[j]+L)2---->转换为一般格式

      b  +      kx        =        y

      对于任意j可以表示为图像上一点J(s[j]+L,f[j]+s[i]2+(s[j]+L)2),目的是求出f[i],因为f[i]的值等于y=kx+b直线的截距。

步骤三:

      问题转化为:对于一条斜率为k=2*s[i]的直线(为什么将它作为斜率,请思考),使其至少过一点J(即上文那个任意J点,j<i),使得该直线的截距最小(其实际意义是为P教授省钱)。那么考虑哪些J点有望成为答案?观赏下图:

                                            image

     首先确定直线斜率的变化特点。由于k=s[i]*2所以我们获得两个重要信息,即k>0并且k随i单调递增。我们将这个斜率画到图上,直线变化如下:

                                            image

     直线斜率的单调递增使得我们能够优化DP。怎么优化?优化的基础是,对于本题,在求最小值的背景下,我们维护一个下凸包,然后随着i的不断循环,凸包的尾巴(横坐标较小部分)会被逐渐削减,同时前端会不断变化。这个过程能够被单调队列完美胜任。去尾巴的原因在下图:

                                       image

·由上文结论可以知道,只要我们的f[i]线过了一个J点,那么表示的是f[j]向 f[i]的状态转移。而我们的目的是使得f[i]线过了一个点并使其的截距最小(当然图中有不合理之处,本题中直线截距应当大于0)。图中可以看出,线段AB的斜率小于f[i]线的斜率,那么意味着f[i]线截距最小值时肯定不会过点A,因为可以通过平移到达更美妙的点B或点C等。所以我们就需要删掉队列中A点,因为后来的所有f[i]斜率必将大于AB,A点将会像凸包内的点一样对答案毫无价值。那么在前端加入点怎么解释呢?那也是同样的道理:

                                    image

·在完成f[i]的状态转移后,i++。那么此时以前的f[i]也就成了f[j]中的一员,我们将它标记为点NEW。这个点是否有资格进入单调队列呢(即是否有资格作为候选的状态转移来源呢)就需要再次利用上文相同的思想。如图,如果 NEW点与A点组成的直线斜率小于线段AB,BC那么我们就要将BC两点从队首弹出,原因是既然有了NEW点的存在,那么以后的直线更愿意选择NEW点,因为这样会使截距更小。

·在代码实现“去除尾巴”和“清理头部”中,都是以判断斜率为标准。其实这类DP代码就是在原来基础上增加了单调队列的两排操作和一个求斜率的函数调用。

·代码之前的一个提醒:上文中的关于j的式子里出现了s[i]2,但不要错误理解为双变量,因为在对于同一个f[i]计算时,求斜率坐标相减就已经直接抵消掉了这个关于i的部分。

 

 1 #include<cstdio>
 2 #define ll long long
 3 #define Empty (head>=tail)
 4 #define go(i,a,b) for(int i=a;i<=b;i++)
 5 const int N=50100;ll s[N],Q[N],f[N],n,x,head,L,tail,j;
 6 inline double X(ll i){return s[i];}
 7 inline double Y(ll i){return f[i]+(s[i]+L-1)*(s[i]+L-1);}
 8 inline double Rate(ll i,ll k){return (Y(k)-Y(i))/(X(k)-X(i));}
 9 int main()
10 {
11     scanf("%lld%lld",&n,&L);s[0]=0;L++;head=1;tail=1;Q[1]=0;
12     go(i,1,n){scanf("%lld",&s[i]);s[i]+=s[i-1];}go(i,1,n)s[i]+=i;
13     go(i,1,n)
14     {
15         while(!Empty&&Rate(Q[head],Q[head+1])<2*s[i])head++;
16         j=Q[head];f[i]=f[j]+(s[i]-s[j]-L)*(s[i]-s[j]-L);
17         while (!Empty&&Rate(Q[tail-1],Q[tail])>Rate(Q[tail],i))tail--;Q[++tail]=i;
18     }
19     printf("%lld\n",f[n]); return 0;
20 }//Paul_Guderian

 

 

[2]仓库建设【LINK】

步骤一:
     列出DP方程式:设f[i]表示在i处建立一个美妙仓库,并且已经安置好i以前的所有货物的最小花费。为了便于计算,进行预处理:
     对于k处货物向i地(i>k)运输,输入时有x[i]表示i到山顶的距离,p[i]表示仓库i货物数量,那么关于运输k处货物花费Wk=(x[i]-x[k])*p[k]。这样看难以进行前缀处理,所以拆开一看,哇:
Wk=x[i]*p[k]-x[k]*p[k]

    我们结合实际情况地定义这些数组来便于操作:令P[i]表示p[1~i]前缀和,令g[i]表示(-x[i]*p[i])的前缀和。我们可以顺利写出Dp方程式:

 f[i]=min(f[j]+x[i]*(P[i-1]-P[j])+g[i-1]-g[j]+c[i]) [0<=j<i]

步骤二:

     对于范围内的j,进行式子变形,获得直线解析式:

    [一句提醒:由上文结论可知,如果变形出来的一部分式子仅与i相关,那么会在计算时被抵消,我们无须在意,但为了便于理解,使用蓝色标记这样的无意义的部分]

     f[i]=f[j]+x[i]*P[i-1]-x[i]*P[j]+g[i-1]-g[j]+c[i]                   

     f[i]=f[j]-x[i]*P[j]-g[j]+x[i]*P[i-1]+g[i-1]+c[i]      

    f[i]+x[i]*P[j]=f[j]-g[j]+x[i]*P[i-1]+g[i-1]+c[i]           

      b +   kx    =              y        

步骤三:

     同上。使用单调队列维护,斜率K=x[i]保持单调递增,满足斜率优化的条件。本题需要额外注意的是f[i]的初始化,它的一个转移来源是把之前所有的货物都收进它的仓库,这个作为f[i]的初值就没问题啦。

代码GOGOGO:

 

 1 #include<cstdio>
 2 #include<algorithm>
 3 #define ll long long
 4 #define X(i) (1.0*p[i])
 5 #define Empty (head>=tail)
 6 #define Y(i) (1.0*f[i]-g[i])
 7 #define go(i,a,b) for(int i=a;i<=b;i++)
 8 const int N=1000006;
 9 int n,tail=0,head=1,Q[N];ll p[N],c[N],f[N],x[N],g[N];
10 double Rate(int i,int j){return (Y(i)-Y(j))/(X(i)-X(j));}
11 int main()
12 {
13     scanf("%d",&n);
14     go(i,1,n)scanf("%d%d%d",x+i,p+i,c+i);
15     go(i,1,n)g[i]=g[i-1]-p[i]*x[i],p[i]+=p[i-1];
16     go(i,1,n)
17     {
18         f[i]=x[i]*(p[i-1])+g[i-1]+c[i];
19         while(!Empty&&Rate(Q[head+1],Q[head])<x[i])head++;
20         int j=Q[head];f[i]=std::min(f[i],f[j]+x[i]*(p[i-1]-p[j])+g[i-1]-g[j]+c[i]);    
21         while(!Empty&&Rate(Q[tail],Q[tail-1])>Rate(i,Q[tail]))tail--;
22         Q[++tail]=i;
23     }
24     printf("%lld\n",f[n]);return 0;
25 }//Paul_Guderian

 

[3]土地购买【LINK】

步骤一:

     列出DP方程式:f[i]=? 咋整啊,这不是序列,方块可以任意选,而且不是连续的!卡机……

步骤O:

      对于这道题,我们想要取得尽量少的总花费,有一个清晰的意识是我们要尽可能利用包含关系。如果一个矩形的长宽比一个矩形都小,那么我们就可以直接忽略!我们的目的求最小的花费,也就是最小的长乘宽的和,可以理解为几个矩形的面积和,如图:

        image

       在去掉被其他矩形包含的无用矩形的情况下,购买这两个矩形“套餐”的价格为:a x b。所以我们发现,蓝色的部分(即相较于原来矩形多出的部分)越小,我们的总共花费就越接近于矩形的本来面积(即图中的黑色和黄色矩形所占据的面积),这意味着我们只要使得矩形的长宽差异越接近,那么多出的蓝色部分花费就越少,由于矩形总面积不变,这样就使得总花费越小。怎样才能使得长宽接近?

     SORT函数解决问题。我们将矩形长作为第一关键字,宽作为第二关键字,全都以从小到大加以排序。这样就使得长宽尽量接近。

     这样做怎样去除包含矩形呢?我们发现对于排序后的矩形i,那么它的长必定大于等于i-1,宽的大小关系不定。如果我们发现i-1的宽也小于矩形i了那么说明这是包含关系——我们就像这样不断向前查找是否具有包含关系,直到i-k宽大于i为止。这样处理后的结果是:序列中的矩形的宽呈现递减关系(否则还会有包含关系),那么这样更加便于我们找到一段区间的宽的最值,DP方程式变得简单。

 步骤一:

             列出DP方程式:在排序后,设f[i]表示处理完前i个矩形的最小花费,且有i处作为一组的结尾。我们使用a[i],b[i]分别表示矩形i的长和宽,根据宽的递减关系,以及a[i]单调递增,方程式如下,表示j是前一组的结尾,j+1到i是新的一个套餐:

      f[i]=min(f[j]+a[i]*b[j+1])    [0<=j<i]

步骤二:

     对于范围内的j,进行式子变形,获得直线解析式:

         f[i]=f[j]+a[i]*b[j+1]

         f[i]-a[i]*b[j+1]=f[j]

         b  +     kx     =y     求直线的截距最小值。

步骤三:

     代码直接就来了:

 

 1 #include<cstdio>
 2 #include<algorithm>
 3 #define ll long long
 4 #define Y(i) 1.0*f[i] 
 5 #define X(i) 1.0*S[i+1].b
 6 #define go(i,a,b) for(int i=a;i<=b;i++)
 7 const int N=50003;
 8 struct P{ll a,b;}S[N];
 9 ll f[N];int n,q[N],t,head=1,tail=1;
10 bool cmp(P A,P B){return A.a==B.a?A.b<B.b:A.a<B.a;}
11 double Rate(int i,int j){return (Y(i)-Y(j))/(X(i)-X(j));}
12 int main(){    
13     scanf("%d",&n);
14     go(i,1,n)scanf("%lld%lld",&S[i].a,&S[i].b);std::sort(S+1,S+n+1,cmp);
15     go(i,1,n){while(t&&S[i].b>=S[t].b)t--;S[++t]=S[i];}n=t;
16     go(i,1,n)
17     {
18         while(head<tail&&Rate(q[head],q[head+1])>=-S[i].a)head++;
19         int j=q[head];f[i]=f[j]+S[i].a*S[j+1].b;
20         while(head<tail&&Rate(q[tail-1],q[tail])<Rate(q[tail],i))tail--;q[++tail]=i;
21     }
22     printf("%lld",f[n]);return 0;
23 }//Paul_Guderian

 

 

[4]特别行动队【LINK】

步骤一:

     列出DP方程式:设f[i]表示处理完前i个士兵的最高战斗力总和,并且在i处分出一组(i作为这一组的结尾)。sum[i]表示1~i的初始战斗力和(前缀)。

 

f[i]=max(f[j]+a*(sum[i]-sum[j])2+b*(sum[i]-sum[j])+c)

                        [0<=j<i]

 

步骤二:

     对于范围内的j,进行式子变形,获得直线解析式:

     f[i]=f[j]+a*(sum[i]-sum[j])2+b*sum[i]-b*sum[j]+c----->展开

     f[i]=f[j]+a*sum[i]2+a*sum[j]2-a*sum[i]*sum[j]+b*sum[i]-b*sum[j]+c

     f[i]+a*sum[i]*sum[j]=f[j]+a*sum[j]2-b*sum[j]+b*sum[i]+a*sum[i]2+c

     b  +        kx       =         y      求直线的截距最大值

和上文不同,这里出现了斜率递减和求最大值,相应的上凸包展示一下:

                       image

步骤三:

     代码直接就来了:

 

 1 #include<cstdio>
 2 #define ll long long
 3 #define X(i) 1.0*sum[i]
 4 #define Empty (head>=tail)
 5 #define go(i,a,b) for(int i=a;i<=b;i++)
 6 #define Y(i) (1.0*f[i]+a*sum[i]*sum[i]-b*sum[i])
 7 const int N=1e6+4;ll a,b,c,x,sum[N],f[N];int Q[N],head,tail,n;
 8 double Rate(int i,int j){return (Y(i)-Y(j))/(X(i)-X(j));}
 9 int main()
10 {
11     scanf("%d%lld%lld%lld",&n,&a,&b,&c);
12     go(i,1,n)scanf("%lld",&x),sum[i]=sum[i-1]+x;
13     go(i,1,n)
14     {
15         while(!Empty&&Rate(Q[head+1],Q[head])>2*a*sum[i])head++;
16         int j=Q[head];f[i]=f[j]+a*(sum[i]-sum[j])*(sum[i]-sum[j])+b*(sum[i]-sum[j])+c;
17         while(!Empty&&Rate(i,Q[tail])>Rate(Q[tail],Q[tail-1]))tail--;
18         Q[++tail]=i;
19     }
20     printf("%lld\n",f[n]);return 0;
21 }//Paul_Guderian

 

 

 

[5]防御准备【LINK】

步骤一:

            列出DP方程式:设f[i]表示在i处放置一个大大的防御塔,并且已经处理好之前的布置的最小花费。已知a[i]为在i处件防御塔的花费,由于这里的两点的距离等于i-j,所以DP方程式有所化简,利用等差数列求和,第一项是1,末项是i-j-1,如下:

    f[i]=min(f[j]+(i-j)*(i-j-1)/2+a[i])  [0<=j<i]

步骤二:

     对于范围内的j,进行式子变形,获得直线解析式:

     f[i]=f[j]+(i-j)*(i-j-1)/2+a[i]------------->同乘2并展开

     2*f[i]=2*f[i]+i2-i*j-i-i*j+j2+j+2*a[i]

     2*f[i]+2*i*j=2*f[j]+j2+j+2*a[i]+i2-i

     b     +  kx =          y            求直线截距的最小值。

步骤三:

     代码美妙就来了:

 1 #include<cstdio>
 2 #define ll long long
 3 #define X(i) (1.0*i)
 4 #define Y(i) (2.0*f[i]+1LL*i*i+i)
 5 #define go(i,a,b) for(int i=a;i<=b;i++)
 6 const int N=1000010;
 7 ll f[N],a[N],Q[N],head,tail,n;
 8 double Rate(int i,int j){return (Y(i)-Y(j))/(X(i)-X(j));}
 9 int main()
10 {
11     scanf("%lld",&n);go(i,1,n)scanf("%lld",&a[i]);
12     go(i,1,n)
13     {
14         while(head<tail&&Rate(Q[head+1],Q[head])<2*i)head++;
15         int j=Q[head];f[i]=f[j]+(1LL*i-j)*(i-j-1)/2+a[i];
16         while(head<tail&&(Rate(i,Q[tail])<Rate(Q[tail],Q[tail-1])))tail--;
17         Q[++tail]=i;
18     }
19     printf("%lld\n",f[n]);return 0;
20 }//Paul_Guderian

 

 

 

[6]序列分割【LINK】

步骤一:

            列出DP方程式:NO!先来一个“预分析”!

     我们发现调皮的小H分割序列方式不是有序的,这样计算分数就很麻烦,我们甚至不能进行任何动态规划操作!我们想个法子,让我们的DP能够顺利从左至右进行下去。尝试发现美丽规律:

image

             这是一个切割后的序列片段。对于序列[1~5],如果我们先切3,再切2和4,那么分数计算为:Score1=(a+b)*(c+d)+a*b+c*d

      还是那个切割后的序列片段。对于序列[1~5],如果我们先切2,再切3,再切4,那么分数计算为:Score2=a*(b+c+d)+b*(c+d)+c*d

     发现美妙结论:Score1=Score2=a*b+a*c+a*d+b*c+b*d+c*d

     (也就是分数总是为两两组合乘积的和,与分割顺序无关!)

那么这样我们就可以像以前一样从容地区间DP了!

     列出DP方程式:设f[i][k]表示在i处砍下第k刀,1~i的区间能够得到的最高分数,并且第k-1刀在j处砍下。使用sum[i]表示元素的前缀和。

    f[i][k]=max(f[j][k-1]+(sum[i]-sum[j])*sum[j])

                    [0<=j<i,0<k<=K]

步骤二:

     对于范围内的j,进行式子变形,获得直线解析式:

     为了降低空间复杂度,可以重新定义f[i]=f[i][k],g[i]=f[i][k-1]这样做本质上就是一个滚动数组。

     f[i]=g[j]+(sum[i]-sum[j])*sum[j]------------->展开

     f[i]=g[j]+sum[i]*sum[j]-sum[j]2-------------->移项

            f[i]-sum[i]*sum[j]=g[j]-sum[j]2-------------->一般化

      b +      kx       =      y     求直线截距的最大值

步骤三:

     代码来临:

 1 #include<cstdio>
 2 #define ll long long
 3 #define X(i) (1.0*sum[i])
 4 #define Empty (head>=tail)
 5 #define Y(i) (1.0*g[i]-sum[i]*sum[i])
 6 #define go(i,a,b) for(int i=a;i<=b;i++)
 7 const int N=100008;ll n,K,sum[N],Q[N],head,tail,f[N],g[N];
 8 double Rate(ll i,ll j){if(X(i)==X(j))return 0;return (Y(i)-Y(j))/(X(i)-X(j));}
 9 int main()
10 {
11     scanf("%d%d",&n,&K);
12     go(i,1,n)scanf("%d",&sum[i]),sum[i]+=sum[i-1];
13     go(k,2,K+1){head=tail=0;go(i,k-1,n)
14     {
15         while(!Empty&&Rate(Q[head+1],Q[head])>-sum[i])head++;
16         int j=Q[head];f[i]=g[j]+sum[j]*(sum[i]-sum[j]);
17         while(!Empty&&Rate(Q[tail-1],Q[tail])<Rate(Q[tail],i))tail--;
18         Q[++tail]=i;}go(j,k-1,n)g[j]=f[j];
19     }
20     printf("%lld",f[n]);return 0;
21 }//Paul_Guderian

 

大米飘香的总结:

     斜率优化的DP在对应的题目中表现出色,并且其完美的优化效果难以被其他思想替代。这也同时决定了它具有很大的局限性——只有上文等能够列出直线解析式并具有决策单调性的情况下才会考虑斜率DP。但无论如何,该思想在其对应的题目领域的作用不可忽视。

    文段主要依靠的是式子简单变形和图形来转化问题。网络上有另一种推导方式,即对于f[i]考虑f[j],f[k]谁对答案的贡献更优,由此列出一个关于j,k的不等式,并加以化简,最终得出维护单调队列的元素进出的判断条件。

     Of course,上文的题目可以归为入门题,这里并没有涉及有关二分等较为复杂的问题。另外,这里还留下两道题作为美妙练习题,代码也给出。

 1 #include<cstdio>
 2 #define ll long long
 3 #define go(i,a,b) for(int i=a;i<=b;i++)
 4 using namespace std;const int N=1000020;
 5 ll f[N],S[N],B[N],a[N],b[N],n;int head,tail,Q[N];
 6 double Rate(int j,int k){return (1.0*f[j]-f[k]+B[j]-B[k])/(1.0*S[j]-S[k]);}
 7 int main()
 8 {
 9     scanf("%lld",&n);
10     go(i,1,n)scanf("%lld",&a[i]);
11     go(i,1,n)scanf("%lld",&b[i]),
12     S[i]=S[i-1]+b[i],B[i]=B[i-1]+b[i]*i;
13     go(i,1,n)
14     {
15         while(head<tail&&Rate(Q[head],Q[head+1])<i)head++;int j=Q[head];
16         f[i]=f[j]+(S[i]-S[j])*i-(B[i]-B[j])+a[i];
17         while(head<tail&&Rate(i,Q[tail])<Rate(Q[tail],Q[tail-1])) tail--;
18         Q[++tail]=i;
19     }
20     printf("%lld\n",f[n]);return 0;
21 }//Paul_Guderian
【BZOJ 3437】
 1 #include<cstdio>
 2 #define go(i,a,b) for(int i=a;i<=b;i++)
 3 const int N=3005;
 4 int A(int x){return x*x;}
 5 int sum[N],g[N],f[N],q[N],n,m,h,r;
 6 double k(int j,int k){return (A(sum[k])-A(sum[j])+g[k]-g[j])/(sum[k]-sum[j]);}
 7 int main()
 8 {
 9     scanf("%d%d",&n,&m);
10     go(i,1,n)scanf("%d",&sum[i]),sum[i]+=sum[i-1],g[i]=A(sum[i]);
11     go(j,1,m-1)
12     {
13         h=r=0;q[0]=j;
14         go(i,j+1,n)
15         {
16             while(h<r&&k(q[h],q[h+1])<2*sum[i])h++;
17             f[i]=g[q[h]]+A(sum[i]-sum[q[h]]);
18             while(h<r&&k(q[r-1],q[r])>k(q[r],i))r--;q[++r]=i;
19         }         
20         go(i,1,n)g[i]=f[i];
21     }
22     printf("%d",m*f[n]-A(sum[n]));return 0;
23 }//Paul_Guderian
【BZOJ 4518】

     Maybe大米饼的思绪很乱。不过是认真的啦。祝你美妙!

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 
 
 
冬天快走了,春天快来了,人们显得特有理想,
我不知道未来怎样,
我将去向何方。—汪峰《我应该真实的生活还是去幻想》
posted @ 2017-07-30 16:23  大米饼  阅读(4835)  评论(9编辑  收藏  举报