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

時間  2019-11-10
標籤 一路走下去的斜率優化動態規劃 简体版
原文   原文鏈接

·隨着網上衆多OIer的步子，大米餅便靜靜地作了如下題目。php

 

·首先列出大米餅的碼風(代碼風格)：前端

①for循環被轉化爲Go循環和Ro循環分別表示升序和降序。②對於維護DP的單調隊列，兩個指針經常使用 Head和Tail兩條。③對斜率優化一類題目的座標點的宏定義X(i)Y(i)，便於理解同時使用double Rate函數計算兩點直線斜率。數組

 

[1]玩具裝箱（詳細闡述） 【LINK】網絡

步驟一：ide

     列出DP方程式:設f[i]表示分組完前i件物品的最小花費,爲方便計算，函數

設sum[i]表示是前i件物品的長度和。優化

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

步驟二：3d

      對於範圍內的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大米餅的思緒很亂。不過是認真的啦。祝你美妙！

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 
 
 
冬天快走了，春天快來了，人們顯得特有理想，
我不知道將來怎樣，
我將去向何方。—汪峯《我應該真實的生活仍是去幻想》
相關文章
相關標籤/搜索
每日一句
    每一个你不满意的现在，都有一个你没有努力的曾经。
本站公眾號
   歡迎關注本站公眾號,獲取更多信息
聯繫我們 最近搜索 最新文章 沪ICP备13005482号-10 MyBatis教程 SQL 教程 MySQL教程 Java 教程 Thymeleaf 教程 Hibernate教程 Spring教程 Redis教程