浅谈斜率优化

对于一类形如：$F_i=\min_{L_i≤j≤R_i} \{F_j+val(i, j) \}$ 的动态规划模型，我们有两种优化方法。

  • 当$val(i,j)$的每一项仅与$i,j$中的一个有关时，我们可以使用单调队列进行优化。
  • 当$val(i,j)$中包含$i,j$的乘积项时，我们可以使用斜率优化

我们主要关注第二种情况，即如何维护斜率以实现快速的转移。

Luogu P2365 任务安排

设：$sumT_i=\sum_{j=1}^{i}T_j,\ sumF_i=\sum_{j=1}^{i}F_j$，$F_i$表示把前$i$个任务分成若干批次处理的最小代价。

不难写出：$F_i=\min_{0≤j<i} \{F_j+sumT_i*(sumF_i-sumF_j)+S*(sumF_N-sumF_j) \}$

如果只是要通过本题，那么$O(n^2)$的算法足矣。但是我们还可以继续优化，最终到达$O(n)$的复杂度。

 

不妨设两个决策$0≤j<k<i$且$k$优于$j$（稍后会解释为什么要设$j<k$）。

有：$F_k+sumT_i*(sumF_i-sumF_k)+S*(sumF_N-sumF_k) <= F_j+sumT_i*(sumF_i-sumF_j)+S*(sumF_N-sumF_j)$

化简移项后有：$F_k-F_j<=(sumT_i-S)(sumF_k-sumF_j)$

注意到$sumF$是单调递增的，我们之前设$j<k$，所以这里$sumF_k-sumF_j>0$（除过去不变号）。我们把它除过去，化成一个斜率式。

即：$\frac{F_k-F_j}{sumF_k-sumF_j}<=SumT_i-S$。

 

注意到这个式子形如：$\frac{y_2-y_1}{x_2-x_1}$这样的斜率式。于是把每个决策$j$看成一个坐标为$(sumF_j,F_j)$的点。

那么上面式子的意义就是：若$j<k$且$j,k$的斜率小于等于$sumT_i-S$，则$k$比$j$优。

那要怎么优化呢？我们考虑维护一个决策队列$q$，注意到如果队列里的决策斜率是单调递增的，那么只要队首$q_l,q_{l+1}$的斜率大于$sumT_i-S$，那么队首就是最优决策。于是我们维护一个斜率单调递增的队列（斜率递增，即一个下凸壳）。

  • 对于队首，我们只需检查$q_l,q_{l+1}$的斜率是否小于等于$sumT_i-S$，如果小于等于，说明$q_{l+1}$优于$q_l$，则将$q_l$出队（由于$sumT-S$递增，所以$q_l$不会再对之后的$sumT_i-S$产生贡献）。
  • 对于队尾，每次转移完$i$后，都要往$q$里加入决策点$(sumF_i,F_i)$。注意到$sumF$也是递增的，所以我们只需比较$(q_{r-1},q_r,i)$这三个点是否满足斜率递增，如果不满足，说明$q_{r-1},q_r$的斜率过大，我们将$q_r$出队，继续判断新的队尾是否构成斜率递增关系即可。

 

每个决策都只会进队和出队一次，所以复杂度是$O(n)$的。

 1 #include<cstdio>
 2 using namespace std;
 3 
 4 const double esp = 1e-8;
 5 const int MAXN = 1000010;
 6 
 7 int N, S, F[MAXN];
 8 int sumT[MAXN], sumC[MAXN];
 9 int l, r, q[MAXN];
10 
11 int X(int i) { return F[i]; }
12 
13 int Y(int i) { return sumC[i]; }
14 
15 double Slope(int i, int k) { return (double) (X(i) - X(k)) / (Y(i) - Y(k)); }
16 
17 int main()
18 {
19     scanf("%d %d", &N, &S);
20     for(int i = 1; i <= N; ++i)
21     {
22         scanf("%d %d", &sumT[i], &sumC[i]);
23         sumT[i] += sumT[i - 1];
24         sumC[i] += sumC[i - 1];
25     }
26     q[l = r = 1] = 0;
27     for(int i = 1; i <= N; ++i)
28     {
29         while(l < r && Slope(q[l + 1], q[l]) - (S + sumT[i]) < esp) ++l;
30         int j = q[l];
31         F[i] = F[j] + sumT[i] * (sumC[i] - sumC[j]) + S * (sumC[N] - sumC[j]);
32         while(l < r && Slope(i, q[r]) - Slope(q[r], q[r - 1]) < esp) --r;
33         q[++r] = i;
34     }
35     printf("%d\n", F[N]);
36     return 0;
37 }

 

拓展：

注意到上面的题目有两个特殊性质，即$sumT-S, sumF$单调递增。第一个性质保证了我们只需保留队列中斜率大于$sumT-S$的部分（就是说从队首出队后，对之后的$sumT-S$都没有贡献。即出队后不可能再入队），第二个性质告诉我们每次新加入的点一定在最末尾，只需和队尾比较斜率即可。

  • 如果$sumT-S$不单调，那我们就不能从队首出队任何元素。由于这是一个单调队列，可以进行二分，在$O(logn)$的时间内找到斜率大于当前的$sumT-S$的决策点位置。
  • 如果$sumF$不单调，我们就需要在任何位置插入决策点，可以用平衡树来维护（即维护一个动态凸壳）。

例题1：Luogu P3628 [APIO2010]特别行动队

令$s_i$为前缀和，$f_i$表示前$i$个士兵的最大和，显然有：$f_i=\max_{0≤j<i} \{f_j+calc(i,j)\}$

$calc(i,j)$表示$j+1~i$的士兵分在一起的代价。

按照套路设$0≤j<k<i$，$k$优于$j$。

于是有：$\frac{(f_k+as_k^2-bs_k)-(f_j+as_j^2-bs_j)}{s_k-s_j}≥2as_i$

维护一个单调递减的队列。发现$2as$是递减的，可以在队首直接进行操纵。因为$s$是递增的，所以一定是在队尾插入决策点。

因为满足上面两个性质，直接$O(n)$维护即可。

 1 #include<bits/stdc++.h>
 2 using namespace std;
 3 
 4 typedef long long ll;
 5 const int MAXN = 1000010;
 6 
 7 int n, a, b, c, l, r;
 8 int x[MAXN], q[MAXN];
 9 ll sum[MAXN], f[MAXN];
10 
11 double slope(int k, int j)
12 {
13     return (double)((f[k] + a * sum[k] * sum[k] - b * sum[k]) - (f[j] + a * sum[j] * sum[j] - b * sum[j])) / (sum[k] - sum[j]);
14 }
15 
16 int main()
17 {
18     scanf("%d %d %d %d", &n, &a, &b, &c);
19     for(int i = 1; i <= n; ++i)
20     {
21         scanf("%d", &x[i]);
22         sum[i] = sum[i - 1] + x[i];
23     }
24     f[0] = q[l = r = 1] = 0;
25     for(int i = 1; i <= n; ++i)
26     {
27         while(l < r && slope(q[l], q[l + 1]) >= 2.0 * a * sum[i]) ++l;
28         f[i] = f[q[l]] + a * (sum[i] - sum[q[l]]) * (sum[i] - sum[q[l]]) + b * (sum[i] - sum[q[l]]) + c;
29         while(l < r && slope(q[r - 1], q[r]) <= slope(q[r], i)) --r;
30         q[++r] = i;
31     }
32     printf("%lld\n", f[n]);
33     return 0;
34 }

 例题2：Luogu P4027 [NOI2007]货币兑换

 可以证明：一定存在一种最优方案，每次兑换金券花光所有的钱，每次换钱时花掉所有的金券。

因为如果第$i$天换成金卷，第$j$天换成钱能够盈利。一定是在第$i$天把钱全部换成金卷，第$j$天再全部换成钱最优。

于是设：$f_i$表示第$i$天的最大收益，$A_i, B_i$分别表示第$i$天最多能得到多少金卷A，B。

那么有：$A_i=\frac{f_iR_i}{a_iR_i+b_i},\  B_i=\frac{f_i}{a_iR_i+b_i},\  f_i=\max_{1≤j<i} \{A_j*a_i+B_j*b_i, f_{i-1}\}$

注意到这里没有单调性，所以我们从另一种角度来推导式子。把转移方程中的$max$和$f_{i-1}$忽略，把关于$j$的项作为横纵坐标。

于是得到：$B_j=-\frac{a_i}{b_i}A_j + \frac{f_i}{b_i}$。可以把式子看作$y=kx+b$的形式。

我们相当于是用一条斜率为$-\frac{a_i}{b_i}$的直线去截每个点，显然当截距最大时，$f_i$最大。

数形结合理解一下，我们从上到下平移这条直线，最优决策点其实构成了一个上凸壳（斜率递减）。

对于当前的斜率$-\frac{a_i}{b_i}$，如果某个点左边的斜率大于它，右边的斜率小于它，这个点就是最优决策。

注意到横坐标，斜率都没有单调性，所以我们还要支持在任意位置插入一个点，即动态维护，查询上凸壳。

我们使用$splay$来维护这个凸壳，$splay$里面按照横坐标排序，快速找到当前点的位置，然后往左右比较，维护斜率的单调递减即可。

  1 #include<bits/stdc++.h>
  2 using namespace std;
  3 
  4 const int MAXN = 100010;
  5 const double EPS = 1e-9;
  6 const double inf = 1e18;
  7 
  8 int n, rt, cnt, fa[MAXN], son[MAXN][2];
  9 double A[MAXN], B[MAXN], lk[MAXN], rk[MAXN], dp[MAXN];
 10 //lk, rk分别表示splay中点k左边，右边直线的斜率
 11 
 12 double calc(int a, int b) {if(fabs(A[a] - A[b]) < EPS) return -inf; return (B[a] - B[b]) / (A[a] - A[b]);}
 13 int chk(int x) {return son[fa[x]][1] == x;}
 14 
 15 void rotate(int x, int &goal)
 16 {
 17     int f = fa[x], g = fa[f], k = chk(x);
 18     if(f == goal) goal = x;
 19     else son[g][chk(f)] = x;
 20     fa[f] = x, fa[x] = g, fa[son[x][k ^ 1]] = f;
 21     son[f][k] = son[x][k ^ 1], son[x][k ^ 1] = f;
 22 }
 23 
 24 void splay(int x, int &goal)
 25 {
 26     while(x != goal)
 27     {
 28         if(fa[x] != goal)
 29         {
 30             if(chk(x) ^ chk(fa[x])) rotate(x, goal);
 31             else rotate(fa[x], goal);
 32         }
 33         rotate(x, goal);
 34     }
 35 }
 36 
 37 int Find(int x, double res)
 38 {
 39     if(!x) return 0;
 40     if(lk[x] + EPS >= res && rk[x] - EPS <= res) return x;
 41     else if(lk[x] - EPS < res) return Find(son[x][0], res);
 42     else return Find(son[x][1], res);
 43 }
 44 
 45 void ins(int &x, int f, int now)
 46 {
 47     if(!x) {x = now, fa[x] = f; return;}
 48     if(A[now] - EPS <= A[x]) ins(son[x][0], x, now);
 49     else ins(son[x][1], x, now); 
 50 }
 51 
 52 int pre(int x)
 53 {
 54     int y = son[x][0], now = y;
 55     while(y)
 56     {
 57         if(lk[y] + EPS >= calc(y, x)) now = y, y = son[y][1];
 58         else y = son[y][0];
 59     }
 60     return now;
 61 }
 62 
 63 int nxt(int x)
 64 {
 65     int y = son[x][1], now = y;
 66     while(y)
 67     {
 68         if(rk[y] <= calc(y, x) + EPS) now = y, y = son[y][0];
 69         else y = son[y][1];
 70     }
 71     return now;
 72 }
 73 
 74 void update(int x)
 75 {
 76     splay(x, rt);
 77     if(son[x][0])
 78     {
 79         int l = pre(x);
 80         splay(l, son[x][0]);
 81         son[l][1] = 0;
 82         lk[x] = rk[l] = calc(x, l);
 83     }
 84     else lk[x] = inf;
 85     if(son[x][1])
 86     {
 87         int r = nxt(x);
 88         splay(r, son[x][1]);
 89         son[r][0] = 0;
 90         rk[x] = lk[r] = calc(x, r);
 91     }
 92     else rk[x] = -inf;
 93     if(lk[x] <= rk[x] + EPS)
 94     {
 95         rt = son[x][0], son[rt][1] = son[x][1];
 96         fa[son[x][1]] = rt, fa[rt] = 0;
 97         lk[rt] = rk[son[rt][1]] = calc(rt, son[rt][1]);
 98     }
 99 }
100 
101 int main()
102 {
103     scanf("%d %lf", &n, &dp[0]);
104     for(int i = 1; i <= n; ++i)
105     {
106         double a, b, r;
107         scanf("%lf %lf %lf", &a, &b, &r);
108         int j = Find(rt, -a / b);
109         dp[i] = max(dp[i - 1], A[j] * a + B[j] * b);
110         B[i] = dp[i] / (a * r + b), A[i] = B[i] * r;
111         ins(rt, 0, i); update(i);
112     }
113     printf("%.3f\n", dp[n]);
114     return 0;
115 }

 

posted @ 2019-09-21 22:13  Aegir  阅读(572)  评论(0)  编辑  收藏  举报