zoukankan      html  css  js  c++  java
  • 斜率优化总结

    好久不见,长假过后焕然一新的我回来了……
    第一次学斜率优化是2年前…现在才写总结【微笑】
    参考:论文《1D1D动态规划优化初步》

    斜率优化是干啥的?

    对动态规划的优化,将其 (O(n^2)) 的复杂度降为 (O(nlogn))(O(n))
    第一次学还是初一,啥都不懂,只是觉得这种类似数形结合的方法真的很妙!

    斜率优化模型

    状态转移方程类似这样:
    $ f[i]=mathop{min}limits_{j=1}^{i-1} lbrace a[i] imes x[j] + b[i] imes y[j] brace $
    我们可以把 (x[i]) 当成横坐标, (y[i]) 当成纵坐标,那么 ((x[i],y[i])) 可当做平面直角坐标系中一个点
    $ f[i]=a[i] imes x[j] + b[i] imes y[j]$ 便表示平面内的一条直线,换个形式写为:
    $ y[j]=-frac{a[i]}{b[i]} imes x[j] + frac{1}{b[i]} imes f[i]$
    其中斜率 (-frac{a[i]}{b[i]}) 为定值,(f[i]) 前系数 (frac{1}{b[i]}) 为定值,对于不同的 (j)(x[j],y[j]) 均已知
    即对每个 (j) 可确定一条直线
    (f[i]) 最小即要这条直线的纵截距最小
    可以理解为一条斜率确定的线从下向上平移,碰到的第一个点为最优决策点 (j)
    而一个重要的性质叫做:所有最优决策点都在平面点集的凸包上
    (为什么大家都说显然啊…画图可直观理解,但证明大概要用线性规划知识…?orz)
    然后可根据这个事实进行优化。

    另一种理解方式

    个人认为可更好想明白为啥在凸包上……
    状态转移方程仍类似这样:
    $ f[i]=mathop{min}limits_{j=1}^{i-1} lbrace a[i] imes x[j] + b[i] imes y[j] brace $
    我们考虑两个决策点 (j)(k) ,若满足 (j>k)(j) 优于 (k),则

    [a[i] imes x[j] + b[i] imes y[j] < a[i] imes x[k] + b[i] imes y[k] \ a[i] imes (x[j]-x[k])<-b[i] imes (y[j]-y[k]) ]

    移项可知 (-frac{a[i]}{b[i]})(frac{y[j]-y[k]}{x[j]-x[k]}) 的大小关系
    (-frac{a[i]}{b[i]}) 是已知确定的,设其为 (kk)
    同样把 ((x[i],y[i])) 看做平面直角坐标系中的点,(frac{y[j]-y[k]}{x[j]-x[k]}) 就为两点所在直线斜率
    上述式子有两种情况:
    情况一:(frac{y[j]-y[k]}{x[j]-x[k]} > kk)
    翻译成文字:若这两点斜率大于 (kk)(j) 优于 (k),反之亦然
    那么所有的最优决策点不会出现如下图情况:

    原因是:图中 (k_{ab}<k_{bc})
    若 $ k_{ab} leq kk$,则 (a) 优于 (b)(b) 不是最优决策点
    (kk<k_{ab}<k_{ac}),则考虑 (bc) 段, (c) 优于 (b)(b) 又不是最优决策点
    所以最优决策点只能在一个斜率递减的凸包上。

    情况二: (frac{y[j]-y[k]}{x[j]-x[k]} > kk)
    与第一种情况类似,(略过长串过程),最优决策点只能在一个斜率递增的凸包上。

    这就是一个很妙的性质啦!

    应用分类

    决策直线斜率与二元组坐标同时满足单调性

    这应该是最经典最常见的应用了。
    由于斜率变化单调,所以决策点在凸包上只会单调移动。
    我们只需维护一个单调队列及决策指针,每次决策时移动指针至最优决策点,决策,然后将新状态插入队列中,删点维护凸壳。
    复杂度 (O(n))
    例题挺多啦,如 (bzoj1010) 玩具装箱,(bzoj1096) 仓库建设……

    这里以 (bzoj3156) 防御准备 为例
    题目:https://www.lydsy.com/JudgeOnline/problem.php?id=3156

    (dp[i]) 表示在第 (i) 个检查点放置守卫塔且前 (i) 个检查点通过安全检查的最小花费
    可列出状态转移方程:

    [egin{equation*} egin{aligned} dp[i]&=mathop{min}limits_{j=0}^{i-1} lbrace dp[j]+sumlimits_{k=j+1}^{i-1} i-k brace+a[i] \ &=mathop{min}limits_{j=0}^{i-1} lbrace dp[j]+frac{(i-j)(i-j-1)}{2} brace+a[i] \ &=mathop{min}limits_{j=0}^{i-1} lbrace dp[j]+frac{1}{2}j^2 + frac{1}{2}j -ij brace+a[i]+frac{1}{2}i^2-frac{1}{2}i \ 其中,dp[0]=0 end{aligned} end{equation*} ]

    (dp[i]+frac{1}{2}i^2 + frac{1}{2}i)(x[i])(i)(y[i])
    原方程可化为 $dp[i]==mathop{min}limits_{j=0}^{i-1} lbrace 1 imes x[j] - i imes y[j] brace+a[i]+frac{1}{2}i^2-frac{1}{2}i $
    最优决策点在斜率递减的凸包上
    决策指针移动时不断删头,最优决策取队首

    代码:
    注意 (dp[0]) 的赋值以及 (long long)

    #include<cstdio>
    #include<iostream>
    #include<algorithm>
      
    using namespace std;
      
    const int N = 1000005;
    typedef long long ll;
      
    ll n;
    ll q[N],h,t;
    ll a[N],dp[N],g[N];
      
    int main()
    {
        scanf("%d",&n);
        for(ll i=1;i<=n;i++) scanf("%lld",&a[i]);
          
        h=t=0;
        dp[0]=g[0]=0;
        q[t++]=0;
        for(ll i=1;i<=n;i++){
            while(t-h>1 && g[q[h+1]]-g[q[h]]<=i*(q[h+1]-q[h])) h++;
            dp[i]=a[i]+1ll*i*(i-1)/2+g[q[h]]-i*q[h];
            g[i]=dp[i]+i*(i+1)/2;
            while(t-h>1 && (g[q[t-1]]-g[q[t-2]])*(i-q[t-1])>=(g[i]-g[q[t-1]])*(q[t-1]-q[t-2])) t--;
            q[t++]=i;
        }
        printf("%lld
    ",dp[n]);
          
        return 0;
    }
    

    其它情况

    一种方法是用平衡树动态维护凸壳
    另一种则是用 (CDQ) 分治。
    复杂度 (O(nlogn))

    例题如 (bzoj1492[NOI2007]) 货币兑换
    题目:https://www.lydsy.com/JudgeOnline/problem.php?id=1492

    题目真是复杂666
    (dp[i]) 表示第 (i) 天可获得的最多的钱
    状态转移方程:

    [egin{equation*} dp[i]=max lbrace dp[i-1],mathop{max}limits_{j=1}^{i-1} lbrace frac{dp[j]}{A[j] imes Rate[j]+B[j]} imes Rate[j] imes A[i] + frac{dp[j]}{A[j] imes Rate[j]+B[j]} imes B[i] brace brace \ 设x[i]=frac{dp[i]}{A[i] imes Rate[i]+B[i]} imes Rate[i],y[i]=frac{dp[i]}{A[i] imes Rate[i]+B[i]} \ 则 dp[i]=max lbrace dp[i-1],mathop{max}limits_{j=1}^{i-1} lbrace x[j] imes A[i] + y[j] imes B[i] brace brace \ end{equation*} ]

    需要维护斜率递减的凸包

    代码:
    先是 (splay) 版的
    细节很多,要多多多注意
    (eps) 要坑死我了。。。

    #include<cstdio>
    #include<iostream>
    #include<algorithm>
    #include<cmath>
     
    #define INF 1e12
    #define eps 1e-9
     
    using namespace std;
     
    const int N = 100005;
    typedef double db;
     
    int n;
    db s,A[N],B[N],R[N],dp[N];
     
    struct node{
        db x,y,lk,rk;
        node *ch[2],*pa;
    }pool[N],*root,*rf,*tmp;
    int cnt;
     
    db sl(node *p,node *q){
        if(fabs(p->x-q->x)<eps) return INF;
        return (p->y-q->y)/(p->x-q->x);
    }
    void rotate(node *p,int ty){
        node *pa=p->pa,*son=p->ch[ty^1],*gp=pa->pa;
        pa->ch[ty]=son; if(son) son->pa=pa;
        p->ch[ty^1]=pa; pa->pa=p;
        p->pa=gp; gp->ch[pa==gp->ch[1]]=p;
        if(root==pa) root=p;
    }
    void splay(node *p,node *q){
        while(p->pa!=q){
            if(p->pa->pa==q)
                rotate(p,p==p->pa->ch[1]);
            else{
                node *pa=p->pa,*gp=pa->pa;
                int f=pa==gp->ch[0]; /**/
                if(p==pa->ch[f]) rotate(p,f),rotate(p,!f);
                else rotate(pa,!f),rotate(p,!f);
            }
        }
    }
    void insert(node *p,node *nd){
        int f=(nd->x>p->x || fabs(nd->x-p->x)<eps);
        if(!p->ch[f]){
            p->ch[f]=nd;
            nd->pa=p;
        }
        else insert(p->ch[f],nd);
    }
    node *findl(node *p,node *nd){ /**/
        if(!p) return p;
        db k=sl(p,nd);
        node *q;
        if(k<p->lk || fabs(k-p->lk)<eps){
            q=findl(p->ch[1],nd);
            return q ? q : p ;
        }
        else return findl(p->ch[0],nd);
    }
    node *findr(node *p,node *nd){
        if(!p) return p;
        db k=sl(p,nd);
        node *q;
        if(k>p->rk || fabs(k-p->rk)<eps){
            q=findr(p->ch[0],nd);
            return q ? q : p ;
        }
        else return findr(p->ch[1],nd);
    }
    void Ins(node *nd){
        insert(root,nd);
        splay(nd,rf);
        node *p=findl(nd->ch[0],nd),*q=findr(nd->ch[1],nd);
        if(p){
            splay(p,root);
            p->ch[1]=NULL;
            nd->lk=p->rk=sl(p,nd);
        }
        else nd->lk=INF;
        if(q){
            splay(q,root);
            q->ch[0]=NULL;
            nd->rk=q->lk=sl(q,nd);
        }
        else nd->rk=-INF;
        if(nd->lk < nd->rk) {// del nd
            if(p && q) { splay(p,rf); splay(q,root); q->ch[0]=NULL; } 
            else if(p) { splay(p,rf); p->ch[1]=NULL; }
            else if(q) { splay(q,rf); q->ch[0]=NULL; }
        }
    }
     
    node *query(node *p,db k){
        if(p->lk<k) return query(p->ch[0],k);
        if(p->rk>k) return query(p->ch[1],k);
        return p;
    }
     
    int main()
    {
        scanf("%d",&n); scanf("%lf",&s);
        for(int i=1;i<=n;i++) scanf("%lf%lf%lf",&A[i],&B[i],&R[i]);
         
        rf=&pool[++cnt];
        for(int i=1;i<=n;i++){
            if(i==1) dp[i]=s;
            else{
                tmp=query(root,-A[i]/B[i]);
                dp[i]=max(dp[i-1],A[i]*tmp->x+B[i]*tmp->y);
            }
            tmp=&pool[++cnt];
            tmp->y=dp[i]/(A[i]*R[i]+B[i]); tmp->x=tmp->y*R[i];
            if(!root) {
                root=tmp;
                root->pa=rf; rf->ch[0]=root;
                root->lk=INF; root->rk=-INF; 
            }
            else Ins(tmp);
        }
        printf("%.3lf
    ",dp[n]);
         
        return 0;
    }
    

    然后是 (CDQ) 分治
    先上一篇 (CDQ) 大神论文,第一个例子就是这题:https://wenku.baidu.com/view/52f9c11cff00bed5b9f31d2d.html

    既然选择了远方,便只顾风雨兼程
  • 相关阅读:
    前端开发中同步和异步的区别
    SQL STUFF函数 拼接字符串
    Javascript的精华
    .net网站报错:对象的当前状态使该操作无效
    免费 WebOffice使用
    DotNetBar教程
    C#获取周的第一天、最后一天、月第一天和最后一天
    C#判断文字是否为汉字
    C# 将字符串转为&#2345;这种的 html实体编码
    SqlServer将没有log文件的数据库文件附加到服务器中
  • 原文地址:https://www.cnblogs.com/lindalee/p/9532183.html
Copyright © 2011-2022 走看看