继四边形不等式和单调队列之后的另一个DP优化方法,这样常见的DP优化方法基本全了。

斜率优化DP

        在DP状态转移方程中,常常可以看到这样的一种形式:

        $\min$取$\max$也是可以的。
        如果$f(i,x)$只与$i$有关或者只与$x$有关直接$O(n)$就能做了。这里麻烦的是转移代价与$x$和$i$都有关系,我们必须遍历所有的情况才能找到最值。如果dp的复杂度为$O(n)$,总的时间复杂度就是$O(n^2)$。
        和其它DP优化类似,斜率优化可以直接降维,将复杂度优化到$O(n)$。看一道入门例题:HDU3507
        状态转移方程很好找。如果规定$S(x)$为前缀和,$dp(x)$为前x个元素的最优解,那么有:

        显然是$O(n^2)$复杂度,必然TLE,考虑斜率优化。对于$j$和$k(j>k)$,如果从$j$转移比从$k$转移更优,那么有:

        整理一下得到:

        由于$S(j)>S(k)$(这里严谨的说是$S(j)\geq S(k)$,不过这个特殊处理之后再提),可以得到:

        令$f(x)=dp(x)+S^2(x)$,则有:

        左边的是一个斜率,斜率优化便基于此。
        每当我们求出一个dp值后,就可以立即求出f的值。这里令$S(x)$为横坐标,$f(x)$为纵坐标,将其放到二维坐标系中,能得到下面的图像(部分):

        很显然,这里有$k(a,b)>k(b,c)$,这能说明什么呢?
        考虑斜率与$2S(x)$的关系,有以下三种关系:

  • $2S(x)>k(a,b)>k(b,c)$,可知b优于a,c优于b。
  • $k(a,b)>2S(x)>k(b,c)$,可知a优于b,c优于b。
  • $k(a,b)>k(b,c)>2S(x)$,可知a优于b,b优于c。

        相等的情况可以合并到上面某一类情况中。观察上面的结果,可以发现此时$b$一定不是最优解。这样的话,$b$点就可以剔除,达到优化的效果。
        由于斜率不能单调减,故我们需要维护一个斜率递增的序列,这样可以保证其中每一个元素都是可以取到最优解的。斜率递增其实就是一个下凸包,如下图所示。

        根据前一篇关于二维凸包的博客,我们可以在线性复杂度下用栈维护凸包。
        可是对于凸包上的这么多点,哪一个才能取到最优解?可以发现,如果上图中$k(a,b)<2S(x)<k(b,c)$,那么$b$点必然取到最优解。因此只需要找到第一个斜率大于$2S(x)$的点即可
        怎么找?可以二分也可以从小到大遍历。当然前者更好,不过后者也不是完全没有用处。注意到本题中$S(x)$递增,那么取到最优解的点也具有单调性,对于之前遍历到的元素,直接剔除即可,因为它们一定不会再参与之后的转移。这一步需要从栈底开始遍历(还要删除栈底元素),显然栈是满足不了需求的,故采用队列去实现。由于斜率递增的缘故,它本质上又是一个单调队列。
        这里填一下上面的坑:求斜率时分母不存在怎么办?由于浮点误差等缘故,通常我们不直接用求斜率的方法,而是将分母乘过去,即判定下面的不等式:

        这样分母为0也没有什么问题了。
        另外注意本题中可能有0这个元素,它会导致对于不同的$i$,$(S(i),f(i))$为同一个点,这对维护凸包会造成影响。解决方法是当斜率相同时,仍然将点出队,也就是说我们维护严格单增的单调队列。示例代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
#include <bits/stdc++.h>

#define N (500000+5)
using namespace std;

int n, m, que[N + 500], f, b;
long long op[N], dp[N];

inline long long get(int x, int y) {
return (dp[x] + op[x] * op[x] - dp[y] - op[y] * op[y]);
}

int main() {
ios::sync_with_stdio(false);
while (cin >> n >> m) {
for (int i = 1; i <= n; i++) {
int x;
cin >> x;
op[i] = x + op[i - 1];//前缀和
}
dp[0] = 0, dp[1] = op[1] * op[1] + m, que[0] = 0, que[1] = 1, f = 0, b = 2;//初始化
for (int i = 2; i <= n; i++) {
while (f + 1 < b) {
if (get(que[f + 1], que[f]) < 2ll * op[i] * (op[que[f + 1]] - op[que[f]]))++f;//删去点
else break;
}
dp[i] = dp[que[f]] + (op[i] - op[que[f]]) * (op[i] - op[que[f]]) + m;//转移
while (f + 1 < b) {
if (1ll * get(i, que[b - 1]) * (op[que[b - 1]] - op[que[b - 2]]) <=//注意是<=
1ll * get(que[b - 1], que[b - 2]) * (op[i] - op[que[b - 1]]))
--b;//维护凸包
else break;
}
que[b++] = i;
}
cout << dp[n] << endl;
}
return 0;
}

示例题目

        继续来看斜率优化DP的例子:APIO2014序列分割
        一开始感觉是道水题,但再一想空间和时间,难度陡升。想要解决本题首先需要看出划分顺序是不影响结果的,这样可以强制从左向右划分,大大简化问题。规定$dp(x,y)$为将前x个数划分成y+1块的最高得分,那么有状态转移方程:

        考虑用斜率优化。假设$j$比$k$更优$(j>k)$,那么有:

        移项,得到:

        令$y(x,y)=dp(x,y)-S^2(x)$,那么就是:

        有$S(j)>S(k)$(当然可能取等,处理方式同上),得到:

        发现不等号方向和上面的入门题反过来,不用想也知道这次是维护上凸包。其它处理方式与入门题基本相同。
        这里有一个优化,根据状态的转移次序,我们可以将dp数组压缩至一维,不过还需要一个数组记第二问答案。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
#include<bits/stdc++.h>

#define GET(x, y) (dp[x]-1ll*op[x]*op[x]-dp[y]+1ll*op[y]*op[y])
#define N 100005
using namespace std;
int n, k, op[N], que[N], f, b, ans[N][205];
long long dp[N], dpt[N];

inline void push(int x) {//入单调队列
while (f + 1 < b) {
if (1ll * GET(x, que[b - 1]) * (op[que[b - 1]] - op[que[b - 2]]) >=
1ll * GET(que[b - 1], que[b - 2]) * (op[x] - op[que[b - 1]]))
--b;
else break;
}
que[b++] = x;
}

inline int front(int x) {
while (f + 1 < b) {
if (1ll * GET(que[f + 1], que[f]) > -1ll * (op[que[f + 1]] - op[que[f]]) * op[x])++f;
else break;
}
return que[f];
}

void solve(int x, int y) {//第二问答案,递归进行
if (y == 0)return;
cout << ans[x][y] << " ", solve(ans[x][y], y - 1);
}

int main() {
scanf("%d%d", &n, &k);
for (int i = 1; i <= n; i++) {
int x;
scanf("%d", &x), op[i] = op[i - 1] + x;
}
for (int ik = 1; ik <= k; ik++) {
f = b = 0, push(ik);
for (int i = ik + 1; i <= n; i++) {
int p = front(i);
ans[i][ik] = p, dpt[i] = dp[p] + 1ll * op[p] * (op[i] - op[p]), push(i);
}
for (int i = ik + 1; i <= n; i++)dp[i] = dpt[i];
}
cout << dp[n] << endl, solve(n, k);
return 0;
}