dp优化略解

· · 算法·理论

#define lg log2下文使用 \lg 代表 \log_2

数据结构优化 dp

数据结构优化 dp,顾名思义,就是用数据结构来优化 dp

单调队列优化

单调队列指的就是元素单调的队列,为了维护这个性质,当插入一个点时就会把队尾不满足单调性质的点全部弹出,然后再把元素放到队尾。

『单调队列主要用于维护两端指针单调不减区间最值』——OI Wiki。

这里举一个例子来帮助理解上面那句话,滑动窗口能用单调队列做的原因就是题目要求的是当前“窗口”的最值,这个"窗口"就是原序列上的一个区间,而且这个区间的左端点和右端点都是单调不减的,因此可以用单调队列优化。

那么单调队列怎么维护两端指针单调不减的区间最值呢?

对于当前要求的区间 [l,r],我们让单调队列中维护的元素就相当于经过如下过程生成的一个单调队列:

如果根据题目要求来设置单调队列的单调性(求区间最小值就维护一个单调递增队列,求区间最大值就维护一个单调递减序列),那么此时队首元素就是答案。因为当我们插入一个元素时只会删除前面和它相比不合法的元素,而这个元素必定是会被插入的,因此当前的答案肯定不会被弹掉,否则它就不是当前的答案了,因为一个元素只能被它后面的元素弹掉。

如果从上一个区间 [l',r'] 转移到 [l,r],满足 l'\le l,r'\le r,那么只需要在上一个区间的单调队列的基础上只需要把 [r'+1,r] 的元素入队,然后再弹掉编号小于 l 的就行了。因为每个元素只会入队和出队一次,所以总的时间复杂度是 \Theta(n) 的。

例题:P3957 跳房子

很明显分数具有单调性,于是直接二分 g,求出此时的最大分数。显然有一种动态规划,设 [l,r] 表示能跳的步数的区间,f_i 表示跳到第 i 个点的最大权值和,则有

f_i=\max_{0\le j<i~\land~x_i-x_j\in[l,r]}\{f_j+s_i\}

最大分数即为

\max_{0\le i\le n}f_i

暴力转移 O(n^2),考虑优化。因为 x_i 单调递增,所以可以转移的 j 所在的区间的两个端点单调不降,因此可以使用单调队列进行优化。具体实现如下:

:::info[代码]

#include <iostream>
#define int long long

using namespace std;

constexpr int N = 5e5 + 10, inf = 0x3f3f3f3f3f3f3f3f;
int n, d, k, a[N], b[N], f[N], q[N], head, tail;

bool check(int mid) {
    int l = max(1ll, d - mid), r = d + mid, res = -inf;
    // [l,r] 表示当前能跳的步数区间
    head = tail = 0, f[0] = 0; // 清空单调队列,初始化
    for (int i = 1, ll = 0; i <= n; i++) {
        // 这里用一个单指针 ll 维护插入
        while (ll < i && a[i] - a[ll] >= l) {
            // 单调队列核心代码,插入元素
            while (head < tail && f[q[tail - 1]] <= f[ll]) tail--;
            q[tail++] = ll;
            ll++;
        }
        while (head < tail && a[q[head]] < a[i] - r) head++;
        if (head == tail) f[i] = -inf;
        else f[i] = f[q[head]] + b[i];
        res = max(res, f[i]);
    }
    return res >= k;
}

signed main() {
    ios::sync_with_stdio(false);
    cin >> n >> d >> k;
    for (int i = 1; i <= n; i++)
        cin >> a[i] >> b[i];
    int l = 0, r = a[n] + 1; // r 取理论极值 + 1
    while (l < r) {
        int mid = (l + r) >> 1;
        if (check(mid)) r = mid;
        else l = mid + 1;
    }
    if (l == a[n] + 1) cout << "-1" << endl;
    else cout << l << endl;
    return 0;
}

:::

单调栈

单调栈类似单调队列,即满足单调性质的栈。插入元素时把不符合单调性质的栈顶弹出,再把元素放到栈顶就行了。

单调队列一般用于寻找左侧/右侧第一个更大/更小元素的问题,当然,它不仅仅局限于这些问题,事实上,和单调队列相比,它更加灵活(题也更毒瘤了),所以我这里通过一道例题来讲解单调栈优化 dp。

例题:CF1407D Discrete Centrifugal Jumps

状态设计十分显然,设 f_i 表示到第 i 座大楼的最少跳跃次数。转移按条件分三类:

第一类转移十分显然,f_i\leftarrow f_{i-1}

第二类转移不好处理,考虑把 \min 拆成两个条件:

对于第二个条件,我们维护一个依次插入了 h_1,\cdots,h_i-1 的由栈顶向栈底严格递增的单调栈,那么所有满足第二个条件的 j 必定在单调栈内,因为一个元素出栈的充要条件就是它的后面(栈顶方向)存在一个大于或等于它的 h_c,那么此时 \displaystyle{\max_{k\in[j+1,i=1]}}\ge h_c\ge h_j,就不满足条件了;反之就满足条件。

对于第一个条件,从栈顶(编号大的一方)的 j 向栈底(编号小的一方)的 j 遍历,则 \displaystyle{\max_{k\in[j+1,i-1]}}h_k 单调不降,所以当有一个时刻 \displaystyle{\max_{k\in[j+1,i-1]}h_k}\ge h_i 时,后面(再往栈底走)的 j 就必定不满足条件了。我们发现上面这个跳出条件和单调栈弹出的跳出条件(就是停止跳出的条件)\displaystyle{h_{top}=\max_{k\in[top,top_0]}h_k>h_i} 是一个包含关系(下面包含上面),所以我们可以在弹栈的时候进行更新,更新次数小于等于弹栈次数。

而第三类转移和第二类的唯一区别就是最大值和最小值和大于小于符号互换了,相对与第二类转移的解法只需要把栈改成由栈顶向栈底严格递减的单调栈就行了。

具体实现看代码,及其简短:

:::info[代码]

#include <iostream>

using namespace std;

constexpr int N = 3e5 + 10;
// sa+sat 是从顶到底单调递减的单调栈,用于维护第三类转移
// sb+sbt 是从顶到底单调递增的单调栈,用于维护第二类转移
int n, a[N], f[N], sa[N], sb[N], sat, sbt;

int main() {
    ios::sync_with_stdio(false);
    cin >> n;
    for (int i = 1; i <= n; i++)
        cin >> a[i];
    f[1] = 0, f[0] = n;
    // 这里注意要把 f[0] 初始化为无穷大
    // 因为本题 f 理论上限为 n-1,所以无穷大取 n
    sa[++sat] = sb[++sbt] = 1;
    for (int i = 2; i <= n; i++) {
        f[i] = f[i - 1] + 1;
        while (sat && a[i] < a[sa[sat]])
            f[i] = min(f[i], f[sa[(sat--) - 1]] + 1);
            // 注意这里是 sat-1,下面也是一样的 sbt-1
            // 因为合法的是这个区间,要从这个区间的左端点的左边转移过去
        while (sat && a[i] == a[sa[sat]]) sat--;
        while (sbt && a[i] > a[sb[sbt]])
            f[i] = min(f[i], f[sb[(sbt--) - 1]] + 1);
        while (sbt && a[i] == a[sb[sbt]]) sbt--;
        sa[++sat] = sb[++sbt] = i;
    }
    cout << f[n] << endl;
    return 0;
}

:::

其它数据结构优化

像线段树、树状数组、K-D Tree、树套树这种数据结构一般都只能根据具体的题目来确定用什么和怎么用,没有什么统一的原则,总之就是涉及到区间更新一个值或一个值更新一个区间时可以用。

例题:P3287 [SCOI2014] 方伯伯的玉米田

首先有一个性质:如果要拔高,那拔高的区间的右端点一定是整排玉米最右边的一个,因为这样不会影响相对大小关系,拔高了还会更优。

根据题意就可以的到动态规划的状态:f_{i,j} 表示考虑前 i 个玉米,已经拔高了 j 次时的最大未拔除的玉米数量,转移方程:

f_{i,j}=\max_{1\le i'<i~\land~0\le j'\le j~\land~0\le h_{i'}+j'\le h_i+j}{f_{i',j'}}+1

下面的条件式子是一个三维偏序的结构,因为这里 f 有两维,所以 CDQ 分治不好用,考虑先用枚举顺序卡掉第一维 i,然后用二维树状数组维护最大值(以你为这里左端点都是从最左边的 10 开始的,所以可以用树状数组)就行了。

实现如下:

:::info[代码]

#include <cstring>
#include <iostream>

using namespace std;

// K 为 a[i]+j 的最大值
constexpr int N = 1e4 + 10, M = 510, K = 5510;
int n, m, a[N], c[M][K], f[N][M], ans;

int query(int x, int y) {
    int res = 0;
    for (int i = x; i; i -= i & -i)
        for (int j = y; j; j -= j & -j)
            res = max(res, c[i][j]);
    return res;
}

void update(int x, int y, int v) {
    for (int i = x; i < M; i += i & -i)
        for (int j = y; j < K; j += j & -j)
            c[i][j] = max(c[i][j], v);
}

int main() {
    ios::sync_with_stdio(false);
    cin >> n >> m;
    for (int i = 1; i <= n; i++) cin >> a[i];
    for (int i = 1; i <= n; i++) {
        for (int j = 0; j <= m; j++)
            f[i][j] = query(j + 1, a[i] + j) + 1;
            // j+1 是为了适应树状数组从 1 开始的要求,下同
        for (int j = 0; j <= m; j++)
            update(j + 1, a[i] + j, f[i][j]), ans = max(ans, f[i][j]);
    }
    cout << ans << endl;
    return 0;
}

:::

CDQ 分治优化

这应该算在数据结构里面吧

CDQ 分治主要是用来处理转移点是一坨的线性动态规划(当然要处理不是线性的也行),这利用了它只用处理一部分对另一部分的贡献的特点,因此可以方便地进行转移,或者利用数据结构辅助转移。

思维很简单,但是应用中的变化很多,灵活性极高。

例题:P5979 [PA 2014] Druzyny

注:为了方便说明,下设

有一个非常显然的暴力动态规划:令 f_i 表示把前 i 个小朋友分组能分的最大组数和对应的方案数,那么有

f_i\leftarrow f_{j}~(i-j\in\text{rng}(j+1,i))

(当 f_i\leftarrow f_j 时,如果 f_i 的最大值等于 f_j 的最大值,那么 f_i 的方案数加上 f_j 的方案数;如果 f_i 的最大值小于 f_j 的最大值,那么令 f_i=f_j

可以清晰地感觉到对于一个 i 来说,j 就是离散的一坨,没有任何规律可循。于是这个时候我们就可以利用 CDQ 分治只用处理一部分到另一部分的贡献的性质,使用 CDQ 分治来解决这道题。

具体来说,假设现在正在计算 f_{[l,m]}f_{[m+1,r]} 的贡献。则 f_i,i\in[m+1,r] 会被 f_{j},j\in[l,m] 更新当且仅当 i-j\in\text{rng}(j+1,i)。为了方便计算,把 \text{rng}(j+1,i) 拆成 \text{rng}(j+1,m)\cap\text{rng}(m+1,i),那么条件可以拆分为

对于条件 1,我们可以在对应区间 [j+\text{mxc}(j+1,m),j+\text{mnd}(j+1,m)] 的左端点和右端点的右边一个分别打上加入和删除操作的标记,到遍历到对应的 i 时执行。对于条件 2,这是一个区间的形式,因此可以使用线段树进行区间查询,至此条件 1 就可以转化为在线段树上的单点修改。

所以这里可以用线段树进行维护,把每个 i 应该做的事情(插入、删除)然后在对应的 i 执行,随后线段树区间查询即可,单次更新复杂度 \Theta(n\lg n)

处理完左对右的贡献,跑一遍 CDQ 分治就行了,总时间复杂度 \Theta(n\lg^2 n)。具体实现看代码:

:::info[代码]

#include <iostream>
#include <vector>
#define endl '\n'

using namespace std;

using pii = pair<int, int>;
using pib = pair<int, bool>;
using state = pii;
// state 用于表示 f 的一个状态:pair<int: 最大组数, int: 方案数>

constexpr int N = 1e6 + 10, mod = 1e9 + 7;
constexpr state basic_state(-1919810, 0);

// 合并两个状态
state operator+(const state& a, const state& b) {
    if (a.first == b.first)
        return state(a.first, (a.second + b.second) % mod);
    return a.first > b.first ? a : b;
}

// 合并两个区间(求交集)
pii operator&(const pii& a, const pii& b) {
    return pii(max(a.first, b.first), min(a.second, b.second));
}

int n;
pii a[N];
state f[N];

struct {
    state e[N << 2];

    void pushup(int u) {
        e[u] = e[u << 1] + e[u << 1 | 1];
    }

    void build(int u = 1, int l = 0, int r = n) {
        e[u] = basic_state;
        if (l == r) return;
        int mid = (l + r) >> 1;
        build(u << 1, l, mid);
        build(u << 1, mid + 1, r);
    }

    void update(int pos, const state& val, int u = 1, int l = 0, int r = n) {
        if (l == r) return void(e[u] = val);
        int mid = (l + r) >> 1;
        if (pos <= mid) update(pos, val, u << 1, l, mid);
        else update(pos, val, u << 1 | 1, mid + 1, r);
        pushup(u);
    }

    state query(int ql, int qr, int u = 1, int l = 0, int r = n) {
        if (ql <= l && r <= qr) return e[u];
        int mid = (l + r) >> 1;
        if (qr <= mid) return query(ql, qr, u << 1, l, mid);
        if (ql > mid) return query(ql, qr, u << 1 | 1, mid + 1, r);
        return query(ql, qr, u << 1, l, mid) + query(ql, qr, u << 1 | 1, mid + 1, r);
    }
} seg; // 单点修改,区间查询的线段树

struct {
    // 一个操作是一个 pair<int:要操作的 j, bool: 插入(1)/删除(0)>
    pib e[N];
    int h[N], ne[N], idx;

    void push(int u, const pib& dat) {
        ++idx;
        e[idx] = dat;
        ne[idx] = h[u];
        h[u] = idx;
    }
} lis; // 使用模拟链表来存每个 i 对应的操作

void cdq(int l, int r) {
    if (l == r) return;
    int mid = (l + r) >> 1;
    cdq(l, mid);
    lis.idx = 0;
    for (int i = mid + 1; i <= r + 1; i++)
        lis.h[i] = 0;
    pii cur(0, n);
    for (int i = mid; i >= l; i--) {
        int ll = max(i + cur.first, mid + 1), rr = min(i + cur.second, r);
        cur = cur & a[i];
        if (ll > rr) continue;
        lis.push(ll, pib(i, true)), lis.push(rr + 1, pib(i, false));
    }
    cur = pii(0, n);
    for (int i = mid + 1; i <= r + 1; i++) {
        for (int j = lis.h[i]; j; j = lis.ne[j]) {
            if (lis.e[j].second)
                seg.update(lis.e[j].first, f[lis.e[j].first]);
            else
                seg.update(lis.e[j].first, basic_state);
        }
        if (i == r + 1) break;
        cur = cur & a[i];
        int ll = max(i - cur.second, l), rr = min(i - cur.first, r);
        if (ll > rr) continue;
        pii t = seg.query(ll, rr);
        if (t.second)
            f[i] = f[i] + pii(t.first + 1, t.second);
    }
    cdq(mid + 1, r);
}

int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);
    cin >> n;
    for (int i = 1; i <= n; i++)
        cin >> a[i].first >> a[i].second;
    f[0] = state(0, 1);
    for (int i = 1; i <= n; i++)
        f[i] = basic_state;
    seg.build();
    cdq(0, n);
    if (f[n].second)
        cout << f[n].first << ' ' << f[n].second << endl;
    else cout << "NIE" << endl;
    return 0;
}

:::

练习题

P2254 [NOI 2005] 瑰丽华尔兹

斜率优化 dp

引入

注:以下用 \text{slope}(P_1,P_2) 表示直线 P_1P_2 的斜率。

让我们先看一道例题:玩具装箱。

简要题意

现在有 n 个玩具,长度分别为 C_i,我们可以把连续的玩具 [l,r] 装进一个容器内,制造容器的代价为

\left(\sum_{k=l}^rC_k+r-l-L\right)^2

其中 L 是常数,n\le5\times 10^4,1\le C_i,L\le10^7

求把所有玩具都装入容器内的最小总代价。

很容易就可以列出动态规划方程:设 f_{i} 表示前 i 个玩具装入任意个容器所需的最小费用,转移方程是

f_i=\min_{0\le j<i}\{f_j+(S_i-S_j+i-j-1-L)^2\}

其中

S_i=\sum_{j=1}^iC_j

边界条件为 f_0=0

这样我们就得到了一个 O(n^2) 的算法,很难通过本题吗,我们需要尝试优化

首先把和 j 无关的项拆出去,为了方便计算,我们设 B_i=S_i+i,L'=L+1,则有

\begin{align} f_i&=\min_{0\le j<i}\{f_j+(B_i-B_j-L')^2\}\\ &=\min_{0\le j<i}\{f_j+B_i^2+B_j^2+L'^2-2B_iB_j-2L'B_i+2L'B_j\}\\ &=\min_{0\le j<i}\{f_j+B_j^2+2(L'-B_i)B_j\}+B_i^2+L'^2-2L'B_i\\ &=\min_{0\le j<i}\{f_j+B_j^2+2(L'-B_i)B_j\}+(B_i-L')^2 \end{align}

我们发现这个转移方程中 B_j 的系数与 i 有关,所以我们不能使用单调队列来优化它。像这种我们就只能使用斜率优化。

斜率优化

对于两个可能的决策点(可能用于更新 f_ij 为当前的决策点,而使上式取到 \min 的就是最优决策点) j_1j_2 满足 1\le j_1<j_2<i,若选择 j_2 要比选择 j_1 更优,那么有

\begin{align} f_{j_1}+B_{j_1}^2+2(L'-B_i)B_{j_1}&>f_{j_2}+B_{j_2}^2+2(L'-B_i)B_{j_2}\\ -2B_iB_{j_1}+f_{j_1}+B_{j_1}^2+2L'B_{j_1}&>-2B_iB_{j_2}+f_{j_2}+B_{j_2}^2+2L'B_{j_2}\\ 2B_i(B_{j_2}-B_{j_1})&>(f_{j_2}+B_{j_2}^2+2L'B_{j_2})-(f_{j_1}+B_{j_1}^2+2L'B_{j_1})\\ \end{align} \begin{align} &\because C_i\ge1\\ &\therefore S_i-S_{i-1}=C_{i}\ge1\\ &\therefore B_i-B_{i-1}=C_i+1\ge2\\ &\therefore\frac{(f_{j_2}+B_{j_2}^2+2L'B_{j_2})-(f_{j_1}+B_{j_1}^2+2L'B_{j_1})}{B_{j_2}-B_{j_1}}<2B_i \end{align}

\begin{align} X_j&=B_j\\ Y_j&=f_j+B_j^2+2L'B_j \end{align}

则上式可以表示为

\frac{Y_{j_2}-Y_{j_1}}{X_{j_2}-X_{j_1}}<2B_i

如果让 P_j=(X_j,Y_j),并把它视为一个点,那么 \frac{Y_{j_2}-Y_{j_1}}{X_{j_2}-X_{j_1}} 就是直线 P_{j_1}P_{j_2} 两点的斜率 \text{slope}(P_{j_1},P_{j_2})j_2 更优的条件就变成了

\text{slope}(P_{j_1},P_{j_2})<2B_i

如果有这样的三个点 A,B,C,对应的 j 分别为 j_1,j_2,j_3

满足 k_1>k_2。设 k_0=2B_i,根据上面推出来的式子,有

分类讨论 k_0k_1,k_2 之间的大小关系,通过上面两条,得

从上面所有情况可以发现,无论如何 j_2 都不可能成为最优决策点,所以可以把 j_2 从候选的最优决策点中移出,只保留 j_1j_3,之后变为这样:

对于每个决策点都可以这样操作,于是最后两两决策点之间的斜率是单调上升的,换句话说,这就形成了一个下凸包

那么我们维护这样一个下凸包有什么用呢?让我们考虑一下最优决策点在这个凸包上的位置。假设这个最优决策点为 j_k,那么和它前面的 j_{k-1}j_{k+1} 相比,它肯定都要更优,根据上面斜率推出来的式子,可以得到

\text{slope}(P_{j_{k-1}},P_{j_k})<2B_i<\text{slope}(P_{j_k},P_{j_{k+1}})

同时因为斜率是单调上升的,那么在 j_k 之前(包括 j_k)的决策点对应的 P 两两相邻点的斜率一定是小于等于 2B_i 的,在 j_k 之后(包括 j_k)的决策点对应的 P 两两相邻点的斜率一定是大于 2B_i 的。这就保证了这个 j_k 是最优决策点。

也就是对一个决策点来说,只要它和它前面的点的斜率小于等于 2B_i,并且它和它后面的点的斜率大于 2B_i,我们就能推出这个决策点一定是最优决策点。

考虑边界情况,如果这个决策点是第一个决策点呢?那么我们就只需要计算它和它后面的决策点的斜率是否大于 2B_i 了;反之如果是最后一个决策点,那么我们就只需要计算它和它前面的斜率是否小于等于 2B_i 了。

到这里,我们已经可以通过维护一个下凸包,并且对于每个 i 在它前面所有 j 组成的下凸包中二分最优决策点,从而实现 O(n\lg n) 的复杂度,在本题 n\le5\times10^4 的数据范围下已经能通过本题。

在这里再补充一些具体实现,因为我们维护的是凸包,所以如果不想在凸包里面再多维护一个 j 的话,我们可以通过 XY 推出 f_i。具体地说,设 j 为最优决策点,有

\begin{align} f_i&=f_j+B_j^2+2L'B_j-2B_iB_j+(B_i-L')^2\\ &=Y_j-2B_iX_j+(B_i-L')^2 \end{align}

另一种理解方式

事实上如果这里令 y=Y_j,x=X_j,k=2B_i,b=f_i-(B_i-L')^2 的话,就变成了

b=y-kx

移一下项,得到

y=kx+b

因为 yx 是随当前决策点 j 的变化而变化的,所以当选定一个决策点时 y=kx+b 就是一条过点 (x,y) 的斜率为 k 的直线(有点抽象),而直线的截距 b 就是我们要最小化的东西。放在坐标系中长这样:

(这是一个任意决策点)

(这是最优决策点)

我们称这条直线为目标直线

当我们选择一个决策点时,事实上就是把目标直线平移直到过决策点对应的二维坐标系上的点,然后它和 y 轴的交点(图中是 I)的纵坐标就是 f_i-(B_i-L')^2,从而实现了最小化。

比较上图中两个决策点,设 P_1 表示一个决策点,P_2 表示 P_1 右边的决策点。

\text{slope}(P_1,P_2)<2B_i,根据代数推出的式子,我们知道这样会让 P_2 代表的决策点优于 P_1,换句话说就是从 P_1 走到 P_2 会让答案(截距)更优,在坐标系中表示为从 P_1 走到 P_2 在目标直线向下的法线方向上移动了距离。图中表示为这样:

其中 \vec u 表示从 P_1P_2 的位移,\vec v 表示目标直线在其法线方向上向下的位移。而 \vec v 又会带来与 y 轴交点的下移,最终使截距 b 最小,得到答案。因此一条线段斜率小于目标直线的斜率会让从它左端点到它右端点时使答案更优。

反之如果斜率小于目标直线斜率,就要从右端点到左端点才能更优了。

于是这样无论 b 是什么,可能的最优决策点就一定在下凸包中。

因此在斜率优化的时候可以把方程变成 b=y-kx 的形式,让 b 表示要最优化的值(一般是 f_i 和一些和 i 相关的常量),让 y 表示只和 j 相关的项,然后让 kx 来表示既有 i 又有 j 的项,变形成 f(i)\times g(j) 的形式。

:::info[\Theta(n\lg n) 的代码]

#include <iostream>
#define int long long

using namespace std;

using pii = pair<int, int>; // 充当点用

double slope(pii a, pii b) { // 计算斜率
    // 因为题目保证了 B_i-B_j>=2,所以不用考虑两个点横坐标相等的情况
    // 在其它题里面可能需要特判这种情况
    // 然后根据纵坐标的大小关系返回 +oo 或 -oo
    return 1. * (b.second - a.second) / (b.first - a.first);
}

constexpr int N = 5e4 + 10;
constexpr double EPS = 1e-6;
pii q[N];
int n, L, b[N], f[N], tail;

signed main() {
    ios::sync_with_stdio(false);
    cin >> n >> L; L++; // 计算 L'
    for (int i = 1; i <= n; i++) {
        cin >> b[i];
        b[i] += b[i - 1];
    }
    for (int i = 1; i <= n; i++) b[i] += i;
    q[tail++] = {0, 0}; // 把基本情况 f_0=0 放进去
    for (int i = 1; i <= n; i++) {
        // 二分最优决策点
        // 最优决策点就是和前面的点的斜率小于等于 2 * b[i] 的最靠右的点
        int l = 0, r = tail - 1;
        while (l < r) {
            int mid = (l + r + 1) >> 1;
            // 注意特判第一个点
            if (mid == 0 || slope(q[mid - 1], q[mid]) - EPS <= 2. * b[i]) l = mid;
            else r = mid - 1;
        }
        // 根据上面说的计算 f[i]
        f[i] = q[l].second - 2 * b[i] * q[l].first + (b[i] - L) * (b[i] - L);
        pii p = {b[i], f[i] + b[i] * b[i] + 2 * L * b[i]};
        // 维护下凸包
        while (tail >= 2 && slope(q[tail - 2], q[tail - 1]) > slope(q[tail - 1], p)) tail--;
        q[tail++] = p;
    }
    cout << f[n] << endl;
    return 0;
}

:::

决策单调性优化

让我们再观察一下题目,我们发现 B_i 是单调增加的,那么 2B_i 也是单调增加的。

对于每次的最优决策点 j 来说,因为 \text{slope}(P_{j-1},P_j)<2B_i,所以随着 2B_i 的单调增加,又因为我们维护的是下凸包,斜率单调上升,所以这意味着最优决策点可能往右移动;因为 \text{slope}(P_j,P_{j+1})>2B_i,所以随着 B_i 的单调增加,又因为下凸包斜率单调上升,所以最优决策点可能需要右移来保证这条性质的成立。

感性理解一下,随着 i 的增加,B_i 单调增加,于是最优决策点 j 也随之要么不动要么往右边移动。这就是所谓的决策单调性,即每次更新 f_i 的最优决策点是在单调往右边移动的(也可能不移动)。

更严谨的证明需要四边形不等式,太过于复杂,这里就不展开了。

总之,在证明了这道题的决策单调性后,可以发现如果一个决策点在当前的最优决策点的左边,那么它永远也不可能成为之后决策中的最优决策点。那么我们根本没有必要维护最优决策点前斜率小于等于 2B_i 的所有决策,可以直接把这一部分全部抛弃掉,并永远也不会使用它。

于是可以用单调队列维护它,我们要保证现在在队列中的点中第一个点就是最优决策点,并且每两个点之间的斜率随着点的编号增加单调增加,即在下凸包中新加进来一个队首为最优决策点的条件。

这样我们就不必进行二分查找了,复杂度从 \Theta(n\lg n) 优化到了 \Theta(n)

放一下 \Theta(n) 的代码(我使用的队列是用 head 指向队列头元素,tail 指向队列末元素的下一个位置,所以初始化时 head=tail=0,可以选择自己认为合适的队列写法):

:::info[\Theta(n) 的代码]

#include <iostream>
#define int long long

using namespace std;

using pii = pair<int, int>; // 充当点用

constexpr int N = 5e4 + 10;
constexpr double EPS = 1e-6;
pii q[N];
int n, L, b[N], f[N], head, tail; // 因为是队列,所以有 head 也有 tail

// 计算斜率,和 O(n lg n) 的一样
double slope(pii a, pii b) {
    return 1. * (b.second - a.second) / (b.first - a.first);
}

signed main() {
    ios::sync_with_stdio(false);
    cin >> n >> L; L++; // 计算 L'
    for (int i = 1; i <= n; i++) {
        cin >> b[i];
        b[i] += b[i - 1];
    }
    for (int i = 1; i <= n; i++) b[i] += i;
    head = tail = f[0] = 0;
    q[tail++] = {0, 0};
    for (int i = 1; i <= n; i++) {
        // 最重要的一行
        // 保证 q[head] 就是最优决策点
        while (tail - head >= 2 && slope(q[head], q[head + 1]) - EPS <= 2 * b[i]) head++;
        f[i] = q[head].second - 2 * b[i] * q[head].first + (b[i] - L) * (b[i] - L);
        pii p = {b[i], f[i] + b[i] * b[i] + 2 * L * b[i]};
        // 维护下凸包
        while (tail - head >= 2 && slope(q[tail - 2], q[tail - 1]) > slope(q[tail - 1], p)) tail--;
        q[tail++] = p;
    }
    cout << f[n] << endl;
    return 0;
}

:::

完结撒花

总结一下,斜率优化可以优化动态规划转移方程中有既有决策点 j 又有当前更新点 i 的问题,然后用二分或者数据结构来维护凸包。如果转移方程还有决策单调性就更好了,可以进一步使用决策单调性优化甚至优化到线性复杂度。

在做斜率优化题的时候有一些需要注意的点:

斜率优化的一些变式

每次插入的点的横坐标单调,目标直线斜率不单调

如果直线斜率 k 不单调,那么就不能保证决策的单调性,就需要维护整个凸包,然后在上面二分答案。

每次插入的点横坐标和目标直线的斜率都不单调

这样直线的斜率就可能出现负数,但是分情况讨论的话最终还是会发现目标直线的 k 和两点间斜率的大小关系和在这两点间转移时的优劣关系。

如果插入的点横坐标不单调的话,单调队列就根本无法维护凸包了,这个时候只能用其它的数据结构来维护,例如说是平衡树。或者也可以用 CDQ 分治离线处理,这个后面会讲。

每次插入的点的横坐标不单调,目标直线斜率单调

和上面一种情况一样,也可以用平衡树维护,如果使用 CDQ 的话还可以少维护一维。

具体来说,和下面的玩具装箱改相比,我们不再需要用二分查询了,直接单调队列就可以(可以做到 n\lg n)。

例题选讲

玩具装箱改

题目链接(完全不保证不保证完全没有锅)

这道题中唯一实质上变了的就是玩具的长度——可以是负的,转移方程并没有改变,于是还按照玩具装箱的处理方式,但这一次 B_i 并不是单调增加的了,因此目标直线的斜率也不是单调增加的了,这符合变式中的第二种情况。

让我们再来分析一下这个 dp 方程

f_i-(B_i-L')^2=\min_{0\le j<i}\{f_j+B_j^2+2(L'-B_i)B_j\}

抽象为

b_i=\min_{0\le j<i}\{y_j-k_ix_j\}

那么对于两个 j_1j_2,若满足 x_{j_1}<x_{j_2},且 j_2j_1 要优,有

\begin{align} y_{j_2}-k_ix_{j_2}&<y_{j_1}-k_ix_{j_1}\\ y_{j_2}-y_{j_1}&<k_i(x_{j_2}-x_{j_1})\\ \frac{y_{j_2}-y_{j_1}}{x_{j_2}-x_{j_1}}&<k_i \end{align}

也就是说,即使可能 j_1>j_2x_{j_1}<x_{j_2},但当我们按 x_j 而不是 j 排序时,斜率优化的结论仍然存在。

(计算斜率的时候注意特判横坐标相等的情况)

那么我们最终仍然是要维护一个下凸包,因为决策单调性的消失,我们只能使用二分;但更大的问题是我们完全无法使用无头单调队列(单调栈)来维护它了。动态维护凸包一般使用平衡树,但如果只是为了 dp 的话可以使用 CDQ。

平衡树做法

在平衡树中以 x 坐标为下标把凸包上的点存进去。在我们插入一个点的时候,可以二分需要删的点的极限位置,也可以暴力找出它的前驱和后继删除。因为每个点最多被加入和删除一次,所以总时间复杂度是 O(n\lg n) 的。查询最优决策点事实上就是找第一个斜率大于目标函数斜率的线段的左端点,所以可以再开一个平衡树来存线段。这种方法可以直接用 set 实现,就不需要手写了。

这个事实上就是计算几何中标准的计算凸包的方法,但是太复杂了,而且我不会,所以不给代码了。

CDQ分治做法

不难发现每个 f_i 都需要前面的 j 来更新,并且另一个限制是需要有序地插入单调队列,这就像是一个二维偏序,因此可以用 CDQ 分治来做。具体做法和 CDQ 分治优化 1D/1D 动态规划(TATT)差不多,就是先计算左边的 dp 值,然后用左边的更新右边的,再计算右边的 dp 值。

:::warning[浮点数的精度误差] 因为无论如何浮点数都有误差,因此我们很有可能因为数据太大而挂掉,于是我们就可以省去计算斜率的这一步,直接计算两条直线斜率的大小关系。这可以用向量叉积来做,也可以简单地两边同乘变为乘法(但是会爆 long long,得开 __int128)。注意还有一大堆特判就是了。

// Pos={__int128,__int128}
// 比较直线 a (过 a1 和 a2) 的斜率 k1 和直线 b (过 b1 和 b2) 的斜率 k2
// 相等返回 0,k1<k2 返回 -1,k1>k2 返回 1
int comp(const Pos& a1, const Pos& a2, const Pos& b1, const Pos& b2) {
    if (a1.x == a2.x) { // 直线 a 垂直于 x 轴
        if (b1.x == b2.x) { // 直线 b 垂直于 x 轴
            // k1=-oo,k2=+oo
            if (a1.y > a2.y && b1.y < b2.y) return -1;
            // k1=+oo,k2=-oo
            else if (a1.y < a2.y && b1.y > b2.y) return 1;
            // k1=k2=+oo/-oo
            return 0;
        } else { // 直线 b 不垂直于 x 轴
            if (a1.y > a2.y) return -1; // k1=-oo<k2
            return 1; // k1=+oo<k2
        }
    } else if (b1.x == b2.x) { // 直线 a 不垂直于 x 轴而直线 b 垂直
        // 和上面差不多
        if (b1.y > b2.y) return -1;
        return 1;
    }
    // 判断斜率相等
    if (((a2.x - a1.x) * (b2.y - b1.y) == (a2.y - a1.y) * (b2.x - b1.x))) return 0;
    // 注意不等式同乘一个数时要记得观察是否需要变方向
    if (((a2.x - a1.x) * (b2.y - b1.y) > (a2.y - a1.y) * (b2.x - b1.x)) ^ ((a1.x < a2.x) ^ (b1.x < b2.x))) return -1;
    return 1;
}

:::

:::info[完整代码]

#include <algorithm>
#include <iostream>
#define int __int128 // 防爆 long long

char wbuf[(1 << 21) + 1], *p3 = wbuf;

#define flush (fwrite(wbuf, 1, p3 - wbuf, stdout), p3 = wbuf)
#define putchar(__x__) (p3 == wbuf + (1 << 21) ? flush : p3, (*p3++) = (__x__))
#define endl putchar('\n')
#define space putchar(' ')

void write(int x) {
    static int stk[100], top;
    if (!x) return void(putchar('0'));
    if (x < 0) putchar('-'), x = -x;
    top = 0;
    while (x) stk[++top] = x % 10, x /= 10;
    while (top) putchar(stk[top--] + '0');
}

void write(const char* str) {
    for (int i = 0; str[i]; i++) putchar(str[i]);
}

char buf[1 << 21], *p1, *p2;
#define getchar() (p1 == p2 && (p2 = (p1 = buf) + fread(buf, 1, 1 << 21, stdin), p1 == p2) ? EOF : *p1++)
int read() {
    int x = 0, f = 1;
    char ch = getchar();
    while (ch < '0' || ch > '9') f = (ch == '-' ? -1 : f), ch = getchar();
    while (ch >= '0' && ch <= '9') x = (x << 1) + (x << 3) + (ch ^ 48), ch = getchar();
    return x * f;
}

using namespace std;

constexpr int N = 2e5 + 10;

struct Pos {
    int x, y;
};

// 斜率比较
int comp(const Pos& a1, const Pos& a2, const Pos& b1, const Pos& b2) {
    if (a1.x == a2.x) {
        if (b1.x == b2.x) {
            if (a1.y > a2.y && b1.y < b2.y) return -1;
            else if (a1.y < a2.y && b1.y > b2.y) return 1;
            return 0;
        } else {
            if (a1.y > a2.y) return -1;
            return 1;
        }
    } else if (b1.x == b2.x) {
        if (b1.y > b2.y) return -1;
        return 1;
    }
    if (((a2.x - a1.x) * (b2.y - b1.y) == (a2.y - a1.y) * (b2.x - b1.x))) return 0;
    if (((a2.x - a1.x) * (b2.y - b1.y) > (a2.y - a1.y) * (b2.x - b1.x)) ^ ((a1.x < a2.x) ^ (b1.x < b2.x))) return -1;
    return 1;
}

Pos q[N]; // 无头单调队列,貌似就是单调栈了
int n, L, b[N], f[N], id[N], tail;

void cdq(int l, int r) {
    if (l == r) return;
    int mid = (l + r) >> 1;
    cdq(l, mid); // 先处理左边的
    // 把两边分别按 x 坐标 (B) 排序
    sort(id + l, id + mid + 1, [&](int i, int j) { return b[i] < b[j]; });
    sort(id + mid + 1, id + r + 1, [&](int i, int j) { return b[i] < b[j]; });
    // 把 [l,mid] 全部扔进单调队列
    tail = 0;
    for (int i = l; i <= mid; i++) {
        Pos p = {b[id[i]], f[id[i]] + b[id[i]] * b[id[i]] + 2 * L * b[id[i]]};
        while (tail >= 2 && comp(q[tail - 1], q[tail], q[tail], p) != -1) tail--;
        q[++tail] = p;
    }
    for (int i = mid + 1; i <= r; i++) {
        // 二分答案
        int ll = 1, rr = tail;
        while (ll < rr) {
            int mid = (ll + rr + 1) >> 1;
            if (mid == 1 || comp(q[mid - 1], q[mid], {0, 0}, {1, 2 * b[id[i]]}) != 1) ll = mid;
            else rr = mid - 1;
        }
        f[id[i]] = min(f[id[i]], q[ll].y - 2 * b[id[i]] * q[ll].x + (b[id[i]] - L) * (b[id[i]] - L));
    }
    // 恢复原状
    for (int i = l; i <= r; i++) id[i] = i;
    // 继续处理右边的
    cdq(mid + 1, r);
}

signed main() {
    n = read(), L = read(); ++L; // 计算 L'
    for (int i = 1; i <= n; i++) b[i] = b[i - 1] + read();
    for (int i = 1; i <= n; i++) b[i] += i;
    for (int i = 0; i <= n; i++) id[i] = i;
    // 答案不会爆 long long,只是运算过程可能会炸
    for (int i = 1; i <= n; i++) f[i] = 0x7fffffffffffffffll;
    cdq(0, n); // 注意 f[0] 会对后面的产生更新
    write(f[n]), endl;
    flush;
    return 0;
}

:::

任务安排

题目链接

让我们先试着列出状态和转移方程。首先很明显设 f_i 表示完成前 i 项任务所需的最小总费用。设

\begin{align} t_i&=\sum_{j=1}^iT_j\\ c_i&=\sum_{j=1}^iC_j \end{align}

但是这里有一个不好解决的问题就是每批任务开始前都需要 s 的时间来初始化机器,在后面的计算中还得加上,这样就会让我们再开一维。但是这 s 的时间带来的额外费用我们是可以提前计算的。具体细节看下面的转移方程:

f_i=\min_{0\le j<i}\{f_j+(c_n-c_j)s+(c_i-c_j)t_i\}

我们把这 s 的启动时间的贡献提前到每次任务中计算。这种优化还是算比较常见的。

很明显这个式子不能单调队列或者单调栈,直接上斜率优化。我们尝试把它化简为 b=y-kx 的形式,满足 b,k 只与 i 有关,y,x 只与 j 有关:

\begin{align} f_i&=f_j+(c_n-c_j)s+(c_i-c_j)t_i\\ f_i&=f_j+c_ns-c_js+c_it_i-t_ic_j\\ f_i-c_ns-c_it_i&=f_j-c_js-t_ic_j \end{align}

\begin{align} y&=f_j-c_js\\ k&=t_i\\ b&=f_i-c_ns-c_it_i\\ x&=c_j \end{align}

则有

\begin{align} b&=y-kx\\ f_i&=y-kx+c_ns+c_it_i \end{align}

观察单调性。因为题目中 T_i 可能为负,所以 t_i,即目标函数斜率,不具有单调性;而 C_i\ge0,所以 c_i 单调不降,即决策点对应的二维平面上的点的横坐标 x 单调不降。符合变式中的第一个,所以维护单调队列,在单调队列上二分即可。

代码不放了,你们应该能写doge。

练习题

斜率优化题单

鸣谢

部分内容来自辰星凌的博客,但改变了一些顺序和描述方式,希望能让你们更好的理解。

文章中可能有错误之处,可以告诉我,望轻喷awa