浅谈序列上的斜率优化 dp

前言

斜率优化是 dp 优化里比较神奇的一派，如果没学过，那么遇到了大概率是不会的，但是只要学过这种思想，那么很多题都可以触类旁通，迎刃而解。

这里主要介绍序列上的斜率优化 dp 模型。

1. 基本形式

dp_i = \min\limits_{L(i)\le j\le R(i)} \text{(或 max)}\{dp_j + val(i,j)\}

其中与单调队列优化不同的是，单调队列优化的针对对象是多项式 val(i,j) 的每一项仅与 i 或 j 有关的情况。

斜率优化针对的是多项式 val(i,j) 包含 i,j 的乘积项的情况。

2. 例题

P2365 任务安排

思路

可见这是个序列问题，考虑将目前处理到第几位作为阶段来划分，即令 dp_{i} 表示把前 i 个任务处理完的最小花费。

枚举最后一次分批在哪个位置，设为 j，0\le j < i。

所以最后一批任务即为 j + 1\sim i。

所以状态转移方程为：

dp_{i} = \min\limits_{0\le j < i}\{dp_{j} + val(j + 1,i)\}

其中 val(j + 1,i) 表示第 j + 1\sim i 个任务带来的花费。

前 j 个任务的最小花费即为 dp_{j}，所以只需计算 val(j + 1, i)。

难点在于花费的计算和完成任务的时间有关，机器启动时间 S 会对后面的阶段造成影响，所以这里需要用到费用提前计算的思想，提前算出启动时间 S 带来的额外花费。

用 sumt_{i} 来表示要完成前 i 个任务所需的时间和，sumc_{i} 表示前 i 个任务的费用系数总和，都可以用前缀和提前求出。

综上所述：状态转移方程为：

dp_{i} = \min\limits_{0\le j < i}\{dp_{j} + sumt_{i}\times (sumc_{i} - sumc_{j}) + S\times (sumc_{n} - sumc_{j})\}

时间复杂度为 O(n^2)，可以通过此题。

\texttt{Code:}

#include <cstring>
#include <iostream>

using namespace std;

const int N = 5010;
typedef long long ll;

int n, S;
ll sumt[N], sumc[N];
ll dp[N];

int main() {
    scanf("%d%d", &n, &S);
    for(int i = 1; i <= n; i++) {
        scanf("%d%d", &sumt[i], &sumc[i]);
        sumt[i] += sumt[i - 1];
        sumc[i] += sumc[i - 1];
    }
    memset(dp, 0x3f, sizeof dp);
    dp[0] = 0;
    for(int i = 1; i <= n; i++)
        for(int j = 0; j < i; j++)
            dp[i] = min(dp[i], dp[j] + sumt[i] * (sumc[i] - sumc[j]) + S * (sumc[n] - sumc[j]));
    printf("%lld", dp[n]);
    return 0;
}

加强版：任务安排2

题目描述一样，但 n\le 3 \times 10^5，考虑优化。

对状态转移方程进行整理，合并同类项，得：

dp_{i} = \min\limits_{0\le j < i}\{dp_{j} - (S + sumt_{i})\times sumc_{j}\} + sumt_{i}\times sumc_{i} + S\times sumc_{n}

对于每个 j，去掉 \min，得：

dp_{i} = dp_{j} - (S + sumt_{i})\times sumc_{j} + sumt_{i}\times sumc_{i} + S\times sumc_{n}

顺便提一句动态规划的优化思路：

一个原则：

在朴素代码上做等价变形。

三个方向：

1. 优化状态设计，阶段和状态转移方程；
2. 若有状态被多次计算或调用，则考虑加上记忆化搜索；
3. 若有多层循环，考虑将与外层循环相关的变量看作定值，及时排除内层循环中的不可能决策或利用数据结构等方法优化找最值的过程，从而优化掉内层循环。

根据上文提到的方向 3，和外层循环 i 相关的值可以暂时看成常量，所以原式的变量就只有 dp_{j} 和 sumc_{j}，分别看成 y 和 x，整理得：

y = (S + sumt_{i})\times x + dp_{i} - sumt_{i}\times sumc_{i} - S\times sumc_{n}

有点眼熟…… 再多看一眼…… 这不是直线斜截式 y = kx + b 吗？

以 sumc_{j} 为横坐标，dp_{j} 为纵坐标，建立平面直角坐标系，所以这是一条斜率为 S + sumt_{i}，截距为 dp_{i} - sumt_{i}\times sumc_{i} - S\times sumc_{n} 的直线。

也就是说，决策候选集合是坐标系的一个点集，每个决策 j 都对应一个点 (sumc_{j}, dp_{j})。

这样一来，状态转移方程就变为了 b_i = \min\limits_{0\le j < i}\{y_i - k_ix_i\}

因为要求的 dp_{i} 是 b_i 中一个系数为正的项，所以问题转化成了：坐标系中有一些点 (x_j,y_j)，从中选出一个点，使得经过它，且斜率为 k_i 的直线的截距最小。

形象一点来说就是有一根斜率为 k_i 的直线从负无穷开始逐渐上移，直到它碰到点集中的点。

如图所示，点 D 就是我们要找的最佳决策点。

观察图像不难发现，只有当 (x_j,y_j) 位于下凸壳时，j 才可能被作为最优决策，所以我们不需要维护整个点集，只需要维护一个下凸壳就行了。

算完 dp_{i} 我们就要把点 (x_i,y_i) 插入点集中，那么如何维护点集才能使我们能快速地找到最优决策呢？

解法一：

首先想到：可以用二分。观察，思考什么时候 (x_j,y_j) 可以作为最优决策呢？

可以发现，当 (x_j,y_j) 与它的上一个点 (x_{j-1},y_{j-1}) 构成的直线的斜率小于 k_i，且与它的下一个点 (x_{j+1},y_{j+1}) 构成的直线的斜率大于 k_i 时，j 为最优决策点。

我们先用一个容器（这里选用队列）来存储目前的点集，然后二分求第一个和上一个点构成的直线斜率小于 k_i 且和后一个点构成的直线斜率大于 k_i 的点。

然后用二分出的答案来计算 dp_{i}。

再插入决策 i，这时要维护下凸壳，检查队尾 q_{tt}、q_{tt - 1} 和要加入的 i 是否构成一个下凸壳，具体地，当

\frac{dp_{q[tt]} - dp_{q[tt - 1]}}{sumc_{q[tt]} - sumc_{q[tt - 1]}} \ge \frac{dp_{i} - dp_{q[tt]}}{sumc_{i} - sumc_{q[tt]}}

时，就应该将队尾弹出。

这样就可以将时间复杂度优化到 O(n\log n)。

\texttt{Code:}

#include <iostream>

using namespace std;

const int N = 300010;
typedef long long ll;
int n, S;
ll sumt[N], sumc[N];
ll dp[N];
int q[N], hh, tt = -1;

int main() {
    scanf("%d%d", &n, &S);
    for(int i = 1; i <= n; i++) {
        scanf("%lld%lld", &sumt[i], &sumc[i]);
        sumt[i] += sumt[i - 1];
        sumc[i] += sumc[i - 1];
    }
    q[++tt] = 0;
    for(int i = 1; i <= n; i++) {
        int l = hh, r = tt;
        while(l < r) {
            int mid = l + r >> 1;
            if(dp[q[mid + 1]] - dp[q[mid]] > (sumt[i] + S) * (sumc[q[mid + 1]] - sumc[q[mid]])) r = mid;
            else l = mid + 1;
        }
        dp[i] = dp[q[l]] - (S + sumt[i]) * sumc[q[l]] + sumt[i] * sumc[i] + S * sumc[n];
        while(hh < tt && (__int128)(dp[q[tt]] - dp[q[tt - 1]]) * (sumc[i] - sumc[q[tt]]) >= (__int128)(dp[i] - dp[q[tt]]) * (sumc[q[tt]] - sumc[q[tt - 1]])) tt--;
        q[++tt] = i;
    }
    printf("%lld", dp[n]);
    return 0;
}

解法二

注意到这道题有个特殊条件：T_i,C_i 都是正整数，意思就是说 sumt 和 sumc 数组一定是单调递增的，所以我们只会不断地在右边加点，而且 k_i 是单调递增的。

外层循环每增加 1，就会新加入一个决策点，同时由于斜率会变大，就可能导致前面的一些决策没用了，要及时排除，只保留凸包上相邻两点连线斜率大于 k_i 的部分，那么凸包最左端的点一定就是最优决策点。

这启示我们可以用单调队列来维护凸包。

建立一个单调队列 q，当外层循环 i 增加 1 后：

1. 排除无用决策，检查队头元素 q_{hh} 和 q_{hh + 1} 构成的直线的斜率是否小于等于当前斜率，具体地，当

\frac{dp_{q[hh + 1]} - dp_{q[hh]}}{sumc_{q[hh + 1]} - sumc_{q[hh]}} \le sumt_{i} + S

时，就应该将队头弹出；

3. 取出队头 j = q_{hh} 为最优决策，用它来计算 dp_{i}；

4. 即将要把 i 插入队尾，检查队尾 q_{tt}、q_{tt - 1} 和要加入的 i 是否构成一个下凸包，具体地，当

\frac{dp_{q[tt]} - dp_{q[tt - 1]}}{sumc_{q[tt]} - sumc_{q[tt - 1]}} \ge \frac{dp_{i} - dp_{q[tt]}}{sumc_{i} - sumc_{q[tt]}}

时，就应该将队尾弹出。

时间复杂度为 O(n)。

\texttt{Code:}

#include <cstring>
#include <iostream>

using namespace std;

const int N = 300010;
typedef long long ll;

int n, S;
ll sumt[N], sumc[N];
int q[N], hh, tt = -1;
ll dp[N];

int main() {
    scanf("%d%d", &n, &S);
    for(int i = 1; i <= n; i++) {
        scanf("%lld%lld", &sumt[i], &sumc[i]);
        sumt[i] += sumt[i - 1];
        sumc[i] += sumc[i - 1];
    }
    q[++tt] = 0;
    for(int i = 1; i <= n; i++) {
        //为防止浮点除法导致精度丢失，故交叉相乘，但要注意极端数据会爆 long long，需使用 __int128
        while(hh < tt && (__int128)(dp[q[hh + 1]] - dp[q[hh]]) <= (__int128)(sumt[i] + S) * (sumc[q[hh + 1]] - sumc[q[hh]])) hh++;
        int j = q[hh];
        dp[i] = dp[j] - sumc[j] * (sumt[i] + S) + sumt[i] * sumc[i] + sumc[n] * S;
        while(hh < tt && (__int128)(dp[q[tt]] - dp[q[tt - 1]]) * (sumc[i] - sumc[q[tt]]) >= (__int128)(dp[i] - dp[q[tt]]) * (sumc[q[tt]] - sumc[q[tt - 1]])) tt--;
        q[++tt] = i;
    }
    printf("%lld", dp[n]);
    return 0;
}

注意：为了防止除法丢失精度，可以交叉相乘。

然鹅，对于一些数据及其毒瘤的题，不仅斜率可能不单调递增，而且横坐标也可能不单调递增，我们分别来看一下这两种情况。

1. 横坐标单增，斜率不单增

P5785 [SDOI2012] 任务安排

题意相同，但 T_i 可能为负数，也就是说 sumt 数组不一定单调递增，也就导致斜率可能一会儿大，一会儿小，单调队列的方法行不通了，就只能采用解法一：二分来做。

2. 横坐标不单增，斜率单增

目前没有遇到这种题。

横坐标不单增，意味着可能需要在中间某个地方插入点，但单调队列还是可以用的，这时可以在单调队列中二分插入，但这样单调队列就要写成“单调链表” 了，感觉不太好写。

3. 横坐标不单增，斜率也不单增

P4655 [CEOI2017] Building Bridges

不得不说这种题真的毒瘤。

这时候就应该拿出斜率优化三件套了：李超线段树、cdq 分治、平衡树。

但是蒟蒻不会前两个，只能选择平衡树。

解释一下为什么是平衡树。

我们需要一个数据结构，支持以下操作：

1. 插入一个点 (x,y)；
2. 查询第一个横坐标小于 x 的点的下标；
3. 查询第一个横坐标大于 x 的点的下标；
4. 删除一个点；
5. 查询第一个和左边点组成的直线斜率小于 k，且和右边点组成的直线斜率大于 k 的点的下标。

综上所述，平衡树是最符合条件的数据结构。

我选择了最灵活的 Splay 来维护凸包。

\texttt{Code:}

#include <cmath>
#include <cstring>
#include <iostream>

using namespace std;

const int N = 100010;
const double eps = 1e-10, inf = 1e15; //注意精度
typedef long long ll;
int n;
ll h[N], sumw[N], dp[N];
struct node{
    int s[2], p;
    double xs, ys, lk, rk;
    void init(int p_, double xs_, double ys_) {p = p_, xs = xs_, ys = ys_;}
}tr[N];
int root, idx;
inline double X(int x) {return h[x];} //横坐标
inline double Y(int x) {return dp[x] + h[x] * h[x] - sumw[x];} //纵坐标
inline double slope(int i, int j) {
    if(fabs(tr[i].xs - tr[j].xs) < eps) return tr[i].ys < tr[j].ys ? inf : -inf; //特判直线垂直于 x 轴的情况
    return (tr[i].ys - tr[j].ys) / (tr[i].xs - tr[j].xs);
}

void rotate(int x) {
    int y = tr[x].p, z = tr[y].p;
    int k = x == tr[y].s[1];
    tr[z].s[y == tr[z].s[1]] = x, tr[x].p = z;
    tr[y].s[k] = tr[x].s[k ^ 1], tr[tr[x].s[k ^ 1]].p = y;
    tr[x].s[k ^ 1] = y, tr[y].p = x;
}

inline void splay(int x, int k) {
    while(tr[x].p != k) {
        int y = tr[x].p, z = tr[y].p;
        if(z != k) {
            if((x == tr[y].s[1]) ^ (y == tr[z].s[1])) rotate(x);
            else rotate(y);
        }
        rotate(x);
    }
    if(!k) root = x;
}

inline void insert(double xs, double ys) {
    int u = root, p = 0;
    while(u) p = u, u = tr[u].s[xs + eps > tr[u].xs];
    u = ++idx;
    if(p) tr[p].s[xs + eps > tr[p].xs] = u;
    tr[u].init(p, xs, ys);
    splay(u, 0);
}

inline int get_prev(int x) {
    int res = 0, u = tr[x].s[0];
    while(u) {
        if (tr[u].lk < slope(u, x) + eps) res = u, u = tr[u].s[1]; //只寻找满足凸包性质的点
        else u = tr[u].s[0];
    }
    return res;
}

inline int get_next(int x) {
    int res = 0, u = tr[x].s[1];
    while(u) {
        if (slope(x, u) < tr[u].rk + eps) res = u, u = tr[u].s[0]; //只寻找满足凸包性质的点
        else u = tr[u].s[1];
    }
    return res;
}
//以上为 Splay 模板

//维护凸壳
inline void update(int &rt, int x) {
    splay(x, 0);
    if (tr[x].s[0]) {
        int y = get_prev(x);
        if(y) { //有可能有左子树但左子树上的点都不满足凸包性质，所以要特判
            splay(y, x);
            tr[y].s[1] = 0; // 清空那些不满足凸包性质的点
            tr[y].rk = tr[x].lk = slope(y, x); //计算斜率
        }
        else tr[x].s[0] = 0, tr[x].lk = -inf; //左边没有点能满足凸包性质，全部删掉且斜率为负无穷
    } else tr[x].lk = -inf;
    if (tr[x].s[1]) { //同上
        int y = get_next(x);
        if(y) {
            splay(y, x);
            tr[y].s[0] = 0;
            tr[x].rk = tr[y].lk = slope(x, y);
        }
        else tr[x].s[1] = 0, tr[x].rk = inf;
    } else tr[x].rk = inf;
    if (tr[x].lk + eps > tr[x].rk) { //若左斜率大于右斜率，则说明不满足凸包性质，删掉
        int y = get_prev(x), z = get_next(x);
        splay(y, 0), splay(z, y);
        tr[z].s[0] = 0;
        tr[y].rk = tr[z].lk = slope(y, z); //重新更新斜率
    }
}

//寻找最优决策
//在 Splay 上二分求解
//若当前斜率在左右斜率之间，则找到
//若当前点的左斜率小于当前斜率，说明要求的点在左侧，递归左子树
//否则在右侧，递归右子树
int get_strategy(int u, double k) {
    if(!u) return 0;
    if(tr[u].lk < k + eps && k < tr[u].rk + eps) return u;
    if(tr[u].lk + eps > k) return get_strategy(tr[u].s[0], k);
    else return get_strategy(tr[u].s[1], k);
}

int main() {
    scanf("%d", &n);
    for(int i = 1; i <= n; i++) scanf("%lld", &h[i]);
    for(int i = 1; i <= n; i++) {
        scanf("%lld", &sumw[i]);
        sumw[i] += sumw[i - 1];
    }
    dp[1] = 0;
    dp[2] = dp[1] + (h[2] - h[1]) * (h[2] - h[1]);
    insert(X(1), Y(1));
    update(root, 1);
    insert(X(2), Y(2));
    update(root, 2); //先插入两个基本决策点，方便计算斜率
    for(int i = 3; i <= n; i++) { 
        int j = get_strategy(root, 2.0 * h[i]); //获取最优决策
        dp[i] = dp[j] + (h[i] - h[j]) * (h[i] - h[j]) + sumw[i - 1] - sumw[j]; //更新
        insert(X(i), Y(i));
        update(root, i); //插入新决策并维护凸包
    }
    printf("%lld", dp[n]);
    return 0;
}

（纯纯毒瘤，调了一个上午 + 下午的一节课时间）

汇总一下各道斜率优化的题：

横坐标单增，斜率也单增

P3195 [HNOI2008] 玩具装箱
P4360 [CEOI2004] 锯木厂选址
P2120 [ZJOI2007] 仓库建设
P3648 [APIO2014] 序列分割
P2900 [USACO08MAR] Land Acquisition G
P4072 [SDOI2016] 征途
P3628 [APIO2010] 特别行动队
AT_dp_z Frog 3

横坐标单增，斜率不单增

P5785 [SDOI2012] 任务安排

横坐标不单增，斜率也不单增

P4655 [CEOI2017] Building Bridges
P4027 [NOI2007] 货币兑换