斜率优化动态规划
推荐论文：
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; }
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; }
