【一路走下去的斜率優化動態規划】


·隨着網上眾多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大米餅的思緒很亂。不過是認真的啦。祝你美妙!

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 
 
 
冬天快走了,春天快來了,人們顯得特有理想,
我不知道未來怎樣,
我將去向何方。—汪峰《我應該真實的生活還是去幻想》


免責聲明!

本站轉載的文章為個人學習借鑒使用,本站對版權不負任何法律責任。如果侵犯了您的隱私權益,請聯系本站郵箱yoyou2525@163.com刪除。



 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM