斜率优化动态规划

推荐论文：

1.从一类单调性看算法的优化汤泽

2.关于DP的斜率优化 http://wenku.baidu.com/view/d3d979dcd15abe23482f4d58.html

3. 题目分析 http://wenku.baidu.com/view/7745777801f69e3143329449.html

4. 动态规划四之四边形不等式和斜率优化 http://wenku.baidu.com/view/383d4de59b89680203d825a3.html

CEOI 2004 锯木厂选址

预处理各前缀和，然后推出普通DP，然后利用决策单调性进行优化，证明过程有点忧伤，不妨跟着论文推推。

#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<algorithm>
using namespace std;

typedef long long LL;
const int inf = 0x3f3f3f3f;
const int N = 20010;
const double esp = 1e-8;
int n;
int d[N], w[N];
int sd[N], sw[N], cost[N];
int F[N];
int Q[N<<1];
int sign(double x){ return x<-esp?-1:(x>esp); }
inline double G(int k1,int k2){
    double A = sw[k1]*sd[k1] - sw[k2]*sd[k2]
    double B = sw[k1] - sw[k2];
    return A/B;
}
inline int ALL(int i,int j){
    return cost[j]-cost[i-1]-sw[i-1]*(sd[j]-sd[i-1]);
}
inline int F(int j,int i){
    return cost[j] + ALL(j+1,i) + ALL(i+1,n+1);
}

int main(){
    while( scanf("%d", &n) != EOF){
        memset( w, 0, sizeof(w));
        memset( d, 0, sizeof(d));
        for(int i = 1; i <= n; i++)
            scanf("%d%d", &w[i], &d[i] );
        w[n+1] = 0; d[n+1] = 0;
        sd[0] = sw[0] = 0;
        for(int i = 1; i <= n+1; i++){
            sd[i] = sd[i-1] + d[i-1];
            sw[i] = sw[i-1] + w[i];
        }
        cost[0] = 0;
        for(int i = 1; i <= n+1; i++){
            cost[i] = cost[i-1] + sw[i-1]*d[i-1];    
        }
        int l = 0, r = 0;
        Q[r++] = 1; Q[r++] = 2;
        f[1] = ALL(2, n+1);            
        f[2] = ALL(3, n+1);    
        for(int i = 3; i <= n; i++){
            while( (l<r) && sign( G(Q[l], Q[l+1]) - sd[i] ) <= 0 )
                ++l;
            f[i] = F( Q[l], i );    
            
            while( (l<r) && (sign( G(Q[r-2],Q[r-1])-G(Q[r-1],i) ) > 0) ) r--;
            Q[r++] = i;
        }
        int ans = inf;
        for(int i = 1; i <= n; i++)
            ans = min( ans, f[i] );
        printf("%d\n", ans );
    }    
    return 0;
}

View Code

hdu 3480 Division

正解应该是四边形优化，不过其决策满足单调性，可利用斜率来优化，2700MS过的.

斜率优化的推理,同上面一题差不多. 算作巩固吧. 在这里写个推理过程:

因为花费是子集中最大减最小的平方, 所以花费最小当然要使最大与最小的差值尽可能小, 所以将 A排序,这个很容易想到.

定义状态 dp( i, j ) 表示前 i个数,将其划分成 j 份,的最小花费. 因为 i+1个数的划分只与前i个数字有关,所以满足无后效性.

所有有转移方程: 其中 j-1 <= k < i

假定对于 (i,j) 而言的最佳决策为k, 有 dp( i, j ) = dp( k, j-1 ) + Ai^2 + A(k+1)^2 - 2*Ai*A(k+1)

令 x = A(k+1) , y = dp( k,j-1) + A(k+1)^2 . 则 dp( i, j ) = y - 2*Ai*x + Ai^2

转换成 dp( i, j ) = Min{ y - 2*Ai*x } + Ai^2;

我们姑且考虑划分集合数量J 为定值, 因为之间没有影响( 只有j-1会影响j, 而相同j之间无影响 ).

我们可以类似上一题假定一个结论: 对于 dp(i,j)的最佳决策为k时, 而dp(i+1,j)的最佳决策为k1, 则必定存在 k1 >= k.

证明如下:

　设 k1 < k2, 而 F( k1 ) = y1 - 2*Ai*x1 , F( k2 ) = y2 - 2*Ai*x2 . // 为了写的方便令 A = 2*Ai

　若 k1 比 k2 更优, 则存在 F( k1 ) < F( k2 ) , 带入计算式有:

　 , 移动下, 因为我们将a排序了,所以有 x1 <= x2 , 所以 x1-x2 <= 0, 除过去括号方向改变,有:

, 令 .

若 g( k1, k2 ) > A , 则有对于 dp( i, j ) 而言, k1 比 k2 更优. 否则 k2 比 k1 优 .

若我们假定决策k 是 dp( i, j )的一个最佳决策, 则存在如下关系:

g( 1..(k-1) , k ) <= A = 2*Ai, 因为 Ai = a(k+1) ,单调非下降.

所以有 g( x, k ) <= 2*Ai <= 2*A(i+1) , 等价于对于 dp( i+1, j ) 而言, 其决策 k 比 (1,k-1) 之间决策都优.

所以我们得到结论, dp( i+1, j ) 的最佳决策 K1 必定与 dp( i, j )的最佳决策 k 存在如下关系, k1 >= k .

证明完毕.

通过以上证明, 我们就可以减少决策数, 来优化DP过程.

其中 j-1 <= k < i

更明白的讲是指, 若我们知道 dp( i, j ) 的最佳决策位于 k, 则 dp( i+1, j ) 的最优决策 k1 则必定是在 [ k, i ] 之间. 就是通过这样来减少决策.起到优化的目的.

到这里. 我们就可以维护一个队列使其满足如下关系:

k1 < k2 < ... < kn .

并且有 g( k1, k2 ) < g( k2, k3 ) < ... < g( kn-1, kn )

维护过程具体如下:

一. 首先利用队中元素更新 dp( i, j ) , 原理如下:

因为队列中保持 k单调, g单调. 上升. 若存在一个 g( k1, k2 ) > 2*Ai. 则说明此时有 F(k1) < F(k2) ,k1比k2更优, 然而 k2以后的 g函数都会大于 2*Ai, 则意味着后面的Kj都不比 k1优.所以 dp( i, j ) 的最佳决策是 k1.

二. 此时将 i位置的信息加入队列, 作为 i+1的一个决策, 此时需要区分 g( kn-1, kn ) 与 g( kn, i ) 的大小关系.

不妨细化来看:

1. g( kn-1, kn ) < g( kn, i ) , 此时若 g( kn-1, kn ) >= 2*Ai. 则对于i+1而言决策优先级有 kn-1 > kn > i . 当前kn-1最优

此时若 g( kn, i ) < 2*Ai , 则对于 I+1而言决策优先级 kn-1 < kn < i ,, 当前I最优

　　　　　　　　　　　　　　　　此时若 g( kn, i ) == 2*Ai, 则优先级: kn-1 > kn , kn < i, 当前最优在 ( kn-1, i )之间.

　　　　　　　　　　　　　　　　　此时若 g( kn, i ) >= 2*Ai, 而 g( kn-1, kn ) < 2*Ai , 则此时有 kn-1 < kn, kn > i . 此时 kn 最优.

　　　　　　　　　　　　　　　　　此情形, 三个决策都可能为最优

2. g( kn-1, kn ) > g( kn, i ), 此时若 g( kn, i ) >= 2* Ai , 则优先级关系: kn-1 > kn > i, 当前 kn-1 最优

此时若 g( kn-1, kn ) < 2*Ai, 则优先级: kn-1 < kn < i, 当前 i 最优

此时若 g( kn-1, kn ) >= 2*Ai, g( kn, i ) < 2*Ai 则最优只可能是 kn-1 或 i.

　　　　　　　　　　　　　　　　此情形下, 最优决策只可能是 kn-1, 或 I.

3. g( kn-1, kn ) = g( kn, i ) , 若 g( kn-1, kn ) >= 2*Ai , 则最优为 kn-1.

　　　　　　　　　　　　　　　　　若 g( kn, i ) < 2*Ai , 则最优为 i.

　　　　　　　　　　　　　　　　　此情形时最优只会是 kn-1, 或者 i.

　　　综上所述 , 当 g( kn-1, kn ) >= g( kn, i ) 时 , 决策 kn 较之 kn-1, i , 不会最优所以可以删去. 否则保留.

整个推理思路就结束了. 然后是代码实现. 考虑下 dp ( i, j ) = dp( k, j-1 ) + Ai ^2 + A(k+1)^2 - 2* Ai * A(k+1)

因为 y = dp( k, j-1 ) + A( k+1 ), x = A(k+1) 我们可以将其转换成 y = dp( k-1, j-1 ) + A(k) , x = A(k) 这样就好写点了.

还有就是当队列不足2个元素时 ,当前只能取一个. 更多细节看看代码吧. 废话有点多,.

#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<algorithm>
using namespace std;

typedef long long LL;
const int N = 10101;
const int M = 5050;
int n, m;

int dp[N][M], a[N], Q[N];

struct Point{
    int x, y;
    Point(){}
    Point(int _x,int _y):x(_x),y(_y){}
}p[N];
int sqr(int x){ return x*x; }

// g(k1,k2) < g(k2,k3)  
int cmp(int a,int b,int c){
    Point p1 = p[a], p2 = p[b], p3 = p[c];
    return (p1.y-p2.y)*(p2.x-p3.x) - (p1.x-p2.x)*(p2.y-p3.y) < 0;
}
// g(k1,k2) >= 2*a_i  
int check(int a,int b,int k){ // a < b
    Point p0 = p[a], p1 = p[b];
    return (p0.y-k*p0.x) <= (p1.y-k*p1.x);  
}
int DP(){
    if( m > n ) return 0;
    for(int i = 1; i <= n; i++)
        dp[i][1] = sqr( a[i]-a[1] );
    for(int j = 2; j <= m; j++){    
        int l = 0, r = 0;
        for(int i = j; i <= n; i++){
            
            p[i] = Point( a[i], dp[i-1][j-1]+sqr(a[i]));            
            while( (l<=r-2) && !cmp(Q[r-2],Q[r-1],i ) ) r--;
            Q[r++] = i;
            while( (l+1<r) && !check(Q[l], Q[l+1], 2*a[i]) ) l++;
            dp[i][j] = p[ Q[l] ].y - 2*a[i]*p[ Q[l] ].x + sqr(a[i]);
        }
    }
    return dp[n][m];
}
int main(){
    int _;
    scanf("%d", &_);
    for(int Case = 1; Case <= _; Case++){
        scanf("%d%d", &n,&m);
        for(int i = 1; i <= n; i++)
            scanf("%d", &a[i] );
        sort( a+1, a+n+1 );    
        printf("Case %d: %d\n", Case, DP() );        
    }
    return 0;
}

View Code