斜率优化dp学习笔记

一年前就看斜率优化dp了…然而一直没有看懂。今天花了一天时间总算了解了个大概。这篇文章将大致分析斜率优化dp的原理和应用。

[HNOI2008]玩具装箱TOY

Description

P教授要去看奥运，但是他舍不下他的玩具，于是他决定把所有的玩具运到北京。他使用自己的压缩器进行压缩，其可以将任意物品变成一堆，再放到一种特殊的一维容器中。P教授有编号为1...N的N件玩具，第i件玩具经过压缩后变成一维长度为Ci.为了方便整理，P教授要求在一个一维容器中的玩具编号是连续的。同时如果一个一维容器中有多个玩具，那么两件玩具之间要加入一个单位长度的填充物，形式地说如果将第i件玩具到第j个玩具放到一个容器中，那么容器的长度将为 x=j−i+∑Ck,i≤K≤j 制作容器的费用与容器的长度有关，根据教授研究，如果容器长度为X,其制作费用为(X−L)2.其中L是一个常量。P教授不关心容器的数目，他可以制作出任意长度的容器，甚至超过L。但他希望费用最小.

Input

• 第一行输入两个整数N，L.接下来N行输入Ci.1≤N≤50000,1≤L,Ci≤107

Output

• 输出最小费用

Sample

Input

5 4
3
4
2
1
4

Output

1

Solution

首先一眼dp方程：

f(i)=min0≤j<n{f(j)+(i−j−1+s[i]−s[j]−L)2}

其中：

s[i]=∑1≤k≤iCi

然而这是个1D/1D方程，如果直接处理复杂度是O(n2)的。状态已经没有办法再简化了，我们只能考虑在转移时做文章。

我们记k(i)为使得f(i)取最小值的j，把方程中对每个i转移的k(i)打一个表，可以发现k(i)是单调不减的。我们考虑证明之。

1. 证明决策单调性

首先对原方程做一些化简。令T(i)=i+s[i],c=L+1，并用f(i,k)表示f(i)在j取k时的值，即：

f(i,j)=f(j)+(T(i)−T(j)−c)2

分析可知要证明：

∀i<j,k(i)≤k(j)

根据k(i)的定义显然有：

∀j,f(i,k(i))≤f(j)

任取i,t满足i<t，只需证：

∀j(j≤k(i)),f(t,k(i))≤f(t,j)

只需证：

f(k(i))+(T(t)−T(k(i))−c)2≤f(j)+(T(t)−T(j)−c)2

设T(t)−T(i)=v，有：

f(k(i))+(T(i)+v−T(k(i))−c)2≤f(j)+(T(i)+v−T(j)−c)2

展开得到：

f(k(i))+(T(t)−T(k(i))−c)2+v2+2v(T(t)−T(k(i))−c)≤f(j)+(T(t)−T(j)−c)2+v2+2v(T(t)−T(j)−c)

由于f(k(i))+(T(t)−T(k(i))−c)2+v2≤f(j)+(T(t)−T(j)−c)2+v2等价于f(i,k(i))≤f(j)，只需证：

2v(T(t)−T(k(i))−c)≤2v(T(t)−T(j)−c)

只需证：

T(k(i))≥T(j)

由于j≤k(i)且T(x)单调增，上式显然成立。

2. 利用单调性解决问题

考虑何时有f(i)取k≥j且取k比取j更优。根据定义：

f(i,k)≤f(i,j)

就是

f(k)+(T(i)−T(k)−c)2≤f(j)+(T(i)−T(j)−c)2

展开并分离变量，得到：

T(i)≥f(k)−f(j)2(T(k)−T(j))+T(k)+T(j)2+c    (!)

这时左边只剩下了和i有关的常量，右边是和k,j有关的代数式。

我们发现这个式子不仅关于k,j对称，而且可以写成y(j)−y(k)x(j)−x(k)+c的形式，前面可以看成斜率的表示。用w(i,j)表示这个斜率，如果a<b<c且w(a,b)>w(b,c)，T(i)≥w(b,c)的条件比T(i)≥w(a,b)更“松”，因而a是不可能被选的。这样“可供选择”的点两两间的斜率形如一个下凸壳，就可以用单调队列维护指标递增-斜率递增的双递增序列。

至此我们成功地获得了一个算法：每次加入新节点时，从队列左端扫去不符合要求的点，剩下的第一个点就是所需的最小答案。之后将当前点插入队列右端并维护下凸性。由于每个元素入队一次出队一次，复杂度为：O(n)。

#include <bits/stdc++.h>
using namespace std;

const int MAXN = 50005;
long long s[MAXN], f[MAXN], L;
int n;
long long q[MAXN];
int l = 1, r = 0;

inline long long T(int i)
{ return i+s[i]; }
long long c;
inline double slop(int j, int k) 
{ return 0.5*(f[k]-f[j])/(T(k)-T(j))+0.5*(T(k)+T(j))+c; }

int main()
{
    scanf("%d%lld", &n, &L);
    c = L+1;
    for (int i = 1; i <= n; i++)
        scanf("%lld", &s[i]), s[i] += s[i-1];
    memset(f, 127/3, sizeof f);
    f[0] = 0;
    q[++r] = 0;
    for (int i = 1; i <= n; i++) {
        while (l < r && slop(q[l], q[l+1]) <= T(i)) l++;
        f[i] = f[q[l]]+(T(i)-T(q[l])-c)*(T(i)-T(q[l])-c);
        while (l < r && slop(q[r-1], q[r]) >= slop(q[r-1], i)) r--;
        q[++r] = i;
    }
    cout << f[n] << endl;
    return 0;
}

四边形不等式

一维情景

在证明决策单调性的过程中计算量非常大，有时可以用四边形不等式简化运算：记j转移到i的代价f(i)−f(j)=w(j→i)。w被称为“凸”的，如果

∀a<b<c<d,w(a→c)+w(b→d)≤w(a→d)+w(b→c)

下面的命题等价：

1. w为凸。
2. w(a→b)+w(a+1→b+1)≤w(a→b+1)+(a+1→b)
3. f决策单调

二维情景

d(i,j)=mini<k<j{d[i,k−1]+d[k+1,j]+w[i,j]}

则下列命题等价：

1. w为凸。
2. w(a→b)+w(a+1→b+1)≤w(a→b+1)+(a+1→b)
3. k(i,j−1)≤k(i,j)≤k(i+1,j)

还记得dp入门题“石子归并”的方程：

d(i,j)=mini≤k<j{d[i,k]+d[k+1,j]+sum[i,j]}

由于sum满足凸性，可以用记忆化搜索记忆k(i,j)将复杂度降到O(n2)。

CEOI2004 锯木厂选址

Description

从山顶上到山底下沿着一条直线种植了n棵老树。当地的政府决定把他们砍下来。为了不浪费任何一棵木材，树被砍倒后要运送到锯木厂。
木材只能按照一个方向运输：朝山下运。山脚下有一个锯木厂。另外两个锯木厂将新修建在山路上。你必须决定在哪里修建两个锯木厂，使得传输的费用总和最小。假定运输每公斤木材每米需要一分钱。

Input

• 输入的第一行为一个正整数n——树的个数（2≤n≤20000）。树从山顶到山脚按照1开始标号。
• 接下来n行，每行有两个正整数（用空格分开）。第i+1行含有：wi——第i棵树的重量（公斤为单位）和 di——第i棵树和第i+1棵树之间的距离，1≤wi≤10000，0≤di≤10000。最后一个数dn，表示第n棵树到山脚的锯木厂的距离。保证所有树运到山脚的锯木厂所需要的费用小于2000000000分。

Output

• 输出只有一行一个数：最小的运输费用。

Sample

Input

9
1 2
2 1
3 3
1 1
3 2
1 6
2 1
1 2
1 1

Output

26 

Solution

“两个”不仅让我们浮想联翩。可以确定一个，枚举另一个作为决策。设dsi为i距离第一棵树的距离，si为前i棵树的总重，定义f(i)：

f(i)=∑1≤k≤jwk(dsj−dsk)+∑j<k≤iwk(di−dk)+∑i<k≤nwk(dn+1−dk)

化简得到：

f(i)=min{Ti+djsj−disj}

其中：Ti=disi+dn+1(sn−si)−∑dkwk。用四边形不等式证明这个函数满足决策单调性：

djsj−disj+dj+1sj+1−di+1sj+1≤djsj−di+1sj+dj+1sj+1−disj+1

化简得到：

di≤di+1

这是显然的。因此可以考虑斜率优化。仿照上面的推导过程可以得到：

di≥dksk−djsjsk−sj

大功告成！原式被整理成了斜率的形式，只需要模仿上面不出脑残错误就可以1A了！(事实上我斜率写错调了1.5h…)

#include <bits/stdc++.h>
using namespace std;

int n;
long long s[20005], d[20005], a[20005], b[20005];

long long q[20005];
int l = 0, r = 0;

double slop(int j, int k) 
{ return 1.0*(d[k]*s[k]-d[j]*s[j])/(s[k]-s[j]); }

int main()
{
    freopen("two.in", "r", stdin);
    freopen("two.out", "w", stdout);
    scanf("%d", &n);
    for (int i = 1; i <= n; i++)
        scanf("%lld%lld", &a[i], &b[i]);
    d[1] = 0;
    for (int i = 2; i <= n+1; i++) d[i] = d[i-1]+b[i-1];
    s[0] = 0;
    for (int i = 1; i <= n; i++) s[i] = s[i-1]+a[i];
    ///////
    long long T = 0;
    for (int i = 1; i <= n; i++) T -= a[i]*d[i]; 
    long long ans = 233333333333ll;
    for (int i = 1; i <= n; i++) {
        while (l < r && slop(q[l], q[l+1]) <= d[i]) l++;
        ans = min(ans, T+d[i]*s[i]-d[n+1]*s[i]+d[n+1]*s[n]-d[i]*s[q[l]]+d[q[l]]*s[q[l]]);
        while (l < r && slop(q[r-1], q[r]) >= slop(q[r-1], i)) r--;
        q[++r] = i;
    }
    cout << ans << endl;
    return 0;
}

参考资料

1. http://www.cnblogs.com/ka200812/archive/2012/08/03/2621345.html
2. http://blog.csdn.net/pi9nc/article/details/9533107
3. http://www.cnblogs.com/kedebug/archive/2013/03/02/2940423.html
4. http://hzwer.com/2114.html
5. https://wenku.baidu.com/view/eeb6d3ea19e8b8f67c1cb937.html
6. 《算法艺术与信息学竞赛（刘汝佳，黄亮）》