0/1 分数规划详解

Jμdge

2018-12-20 16:48:57

Personal

最近入门了 0/1 分数规划,于是写篇博客纪念 然后骗访客量→_→[my cnblogs](https://www.cnblogs.com/Judge/p/10173795.html)(就点一下就好啦!) 分数规划是一项不常用到的(但还蛮实用的)算法,一般来讲就是求一个最优比率。 ## 分数规划的简单介绍 分数规划顾名思义就是求一个分数表达式的最大(小)值。 比如说有 n 个物品,每个物品有两个权值 a 和 b ,然后让你选出**任意件数**(但可能会有限制)的物品,使得两个**权值和**间的比值最大,即求 $\dfrac{\sum_{i=1}^{k} a[i]}{\sum_{j=1}^{k} b[j]}$ (在这里 $1-k$ 为挑出的 k 件物品)的最大值,然后对选择物品方面给出**一定的限制条件**,那么一道 0/1 分数规划的题目就完成了 ## 分数规划的求解方法 分数规划有两种求解方法,一种是 **二分法**,而另一种则是 **Dinkelbach 算法**,一般来讲我们选用第一种方法进行分数规划求解 ### 1.二分法 问题同上,求 $\dfrac{\sum_{i=1}^{k} a[i]}{\sum_{j=1}^{k} b[j]}$ 的最大值,然后加上一个限制:$k>=W$ 我们令 $sum=\sum_{i=1}^{k} a[i] ,\ tot=\sum_{j=1}^{k} b[j]$ ,$k>=W$ 然后假设问题中的最优解为 $ans$ ,那么必然有: ## $$\dfrac{sum}{tot}<=ans$$ 移项得: ## $$sum<=ans\times tot$$ 继续移就得到: ## $$sum-ans\times tot<=0$$ 这样转化有什么用呢?那我们尝试将 $sum$ 和 $tot$ 带回去,就可以得到这么一个式子: ## $$\sum_{i=1}^{k} (a[i]-b[i]\times ans) <=0$$ 这个式子不难理解,就是把整体的贡献转化为了单件物品的贡献。 那么我们只需要二分这个 ans, 计算出每件物品的 $a-b\times ans$,然后排个序,贪心取前 $W$ 个加起来,看看最后的值是否 $<=0$ ,然后就可以根据结果移动左右边界了。 ### 2.Dinkelbach 算法 这个算法其实就是以二分函数图像的单调性为原理,利用了check中计算得出的结果,通过**改变二分枚举的中间点在图像上的位置**来优化二分。 然后先给一张图: ![](https://img-blog.csdn.net/20170206231653349?watermark/2/text/aHR0cDovL2Jsb2cuY3Nkbi5uZXQvaHpvaV96dHg=/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/SouthEast) (途中的红线 x=7 即当前二分的中间点) 在这里我们看到每次枚举的中间点都在一条直线上(但注意单一的直线并不是二分函数图像),那么对于枚举出的中间点其实可以不着急扔,先计算出当前点所在的直线,然后找到这条直线在 x 轴上的截距,用这个截距作为下一次二分的中间点,这样可能会更快逼近正确答案 尽管理解起来不难,但还是配套图吧。下面是算法实现的图组 ![](https://cdn.luogu.com.cn/upload/image_hosting/pgnpwm0t.png) 这里显然红点就是答案 ![](https://cdn.luogu.com.cn/upload/image_hosting/627gmd6b.png) 这里我们第一次可以直接二分中点,得到一个初始点 ![](https://cdn.luogu.com.cn/upload/image_hosting/zkg0o4fd.png) 然后我们过这个得到的点做一条直线,与函数图像相切, 接着我们将该直线与 $x$ 轴交点的 $x$ 坐标作为新点的 $x$ 坐标(好像离答案更远了啊) ![](https://cdn.luogu.com.cn/upload/image_hosting/nlt9juvc.png) 然后新点重复上述过程,做直线 ![](https://cdn.luogu.com.cn/upload/image_hosting/goq4n2ke.png) 这时我们发现已经得到答案了 然鹅这个算法的效果却不见的有多么优秀(尽管上面的例子中这个算法仅仅用了两次),因为二分的函数图像有可能是**坑坑洼洼的**,所以有时这个算法跑数据的时间反而比二分大,但是如果图像是**较为光滑的弧线**,或许$Dinkelbach$ 算法就能充分展现它的优势了。 另外这个算法复杂度就是 $O($玄学$)$ 的,出题人一般不会卡(想卡也难卡,因为初始点可以是任意的,这就相当于随机算法了啊!鬼知道你乱搞出来的初始点会不会就是答案了),但是讲道理,这个算法效率确实没有保障... 至于其他的地方(包括证明),和二分其实差不多,就不啰嗦了 ## 分数规划的有机结合 分数规划一般来讲不会单独成题,一般来讲有以下几种形式: >0.不与任何算法结合,即分数规划裸题 >1.与$0/1$背包结合,即最优比率背包 >2.与生成树结合,即最优比率生成树 >3.与负环判定结合,即最优比率环 >4.与网络流结合,即最大密度子图 >5.与费用流结合,即最优比率流(这个是我乱叫的) >6.与其他的各种带选择的算法乱套,即最优比率啥啥的... 对于上面的后三个结合图论算法,其实类似是把供选择的物品变为了图中的边,然后限制条件求解最优比率。 ### 0.分数规划裸题 这个在介绍二分算法的时候讲过了 #### 题目 [Dropping tests](http://poj.org/problem?id=2976) #### 代码 ```cpp //by Judge #include<cstdio> #include<cstring> #include<iostream> #include<algorithm> #define eps 1e-4 using namespace std; const int M=1005; #ifndef Judge #define getchar() (p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<21,stdin),p1==p2)?EOF:*p1++) #endif char buf[1<<21],*p1=buf,*p2=buf; inline int read(){ int x=0,f=1; char c=getchar(); for(;!isdigit(c);c=getchar()) if(c=='-') f=-1; for(;isdigit(c);c=getchar()) x=x*10+c-'0'; return x*f; } int n,k; double a[M],b[M],c[M],l,r=1e6,mid,sum; inline bool check(){ for(int i=1;i<=n;++i) c[i]=a[i]-b[i]*mid; sort(c+1,c+1+n),sum=0; for(int i=k;i<=n;++i) sum+=c[i]; return sum>=0; } int main(){ while(1){ n=read(),k=read(); if(!n&&!k) return 0; l=0,r=1e6,++k; for(int i=1;i<=n;++i) a[i]=read(); for(int i=1;i<=n;++i) b[i]=read(); while(r-l>eps){ mid=(l+r)/2; check()?l=mid:r=mid; } printf("%.0lf\n",100.0*l); } return 0; } ``` ### 1.最优比率背包 题目大意和裸题的差不多,也是求 n 个物品中,$\dfrac{\sum_{i=1}^{k} a[i]}{\sum_{j=1}^{k} b[j]}$ 的最大值,但是条件有变化,这次我们要对权值 b 进行限制,首先我们还是令 $sum=\sum_{i=1}^{k} a[i] ,\ tot=\sum_{j=1}^{k} b[j]$ ,那么限制条件就是 $tot>=W$($W$已给出) 那么这个时候我们设完 $ans$ 之后得到的式子同样是: ## $$\sum_{i=1}^{k} (a[i]-b[i]\times ans) <=0$$ 于是 $ans$ 照常二分,对于条件成立的判断方式变一变就好了。 其实你应该也猜到了判断方式,原先我们不是贪心取前面的 k 个数么,那这次我们把贪心改成背包就好了啊 那样的话,对于每个物品来说,体积就是 b ,而价值就是上面的 $a[i]-b[i]\times ans$,然后我们跑 0/1 背包,判断满背包价值的正负性。 #### 题目 [Talent Show](https://www.luogu.org/problemnew/show/P4377) #### 代码 ```cpp //by Judge #include<cstdio> #include<cstring> #include<iostream> #define ll long long using namespace std; #ifndef Judge #define getchar() (p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<21,stdin),p1==p2)?EOF:*p1++) #endif inline void cmax(ll& a,ll b){if(a<b)a=b;} inline int Min(int a,int b){return a<b?a:b;} char buf[1<<21],*p1=buf,*p2=buf; inline int read(){ int x=0,f=1; char c=getchar(); for(;!isdigit(c);c=getchar()) if(c=='-') f=-1; for(;isdigit(c);c=getchar()) x=x*10+c-'0'; return x*f; } int n,l,r=1e6,mid,V,t[305]; ll w[305],f[10005]; inline bool check(int tim){ memset(f,-0x3f,sizeof(f)),f[0]=0; ll tmp=f[V]; for(int i=1,j,k;i<=n;++i) for(j=V;~j;--j) if(f[j]!=tmp){ cmax(f[Min(j+w[i],V)],f[j]+t[i]-w[i]*tim); } return f[V]>=0; } inline int calc(){ while(l<=r){ mid=l+r>>1; check(mid)?l=mid+1:r=mid-1; } return r; } int main(){ n=read(),V=read(); for(int i=1;i<=n;++i){ w[i]=read(), t[i]=read()*1000; } return !printf("%d\n",calc()); } ``` ### 2.最优比率生成树 对于构造最优比率生成树的分数规划,其实差不多类似是将 $n$ 物品转换为 $m$ 条边,然后求解,只不过限制改了,选出的 $n-1$ 条边要构成一棵树 具体点说,题目一般会给你一张图(有时候是完全图),然后给你 n 个点 m 条边,求出原图包含 n 个点的一棵树,使得**所选的边的两个权值和比值**最大 那么求法呢,其实也很类似: >>我们考虑构造最小生成树的时候,边是从小到大加的 >>那么对于每一条边,我们让他的边权为 $a[i]-b[i]\times ans$ >>然后和克鲁斯卡尔一样排个序一条一条尝试加就好了(然鹅有时候你得尝试普里姆算法) >>如果加成功了就累计入 $sum$ >>最后根据 $sum$ 的正负性来判断 $check$ 的返回值 其实和前面两个没高明到哪里去,基本想法一样的 #### 题目 [Desert King](http://poj.org/problem?id=2728) #### 代码 (慢着,kruskal会炸?还是我打开方式不对?只好开 $Prim$ ) ``` //by Judge #include<cmath> #include<cstdio> #include<cstring> #include<iostream> #define eps 1e-5 using namespace std; const int M=1005; #ifndef Judge #define getchar() (p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<21,stdin),p1==p2)?EOF:*p1++) #endif char buf[1<<21],*p1=buf,*p2=buf; inline bool cmin(double& a,double b){return a>b?a=b,1:0;} inline int read(){ int x=0,f=1; char c=getchar(); for(;!isdigit(c);c=getchar()) if(c=='-') f=-1; for(;isdigit(c);c=getchar()) x=x*10+c-'0'; return x*f; } int n,pat,x[M],y[M],z[M],vis[M]; double f[M][M],g[M][M],minv[M]; double l,r,mid,sum,tmp; inline bool check(){ memset(vis,0,sizeof(vis)); sum=0,vis[1]=1; for(int i=1;i<=n;++i) minv[i]=g[1][i]-f[1][i]*mid; for(int i=2,j,k;i<=n;++i){ for(tmp=1e16,k=-1,j=2;j<=n;++j) if(!vis[j]&&cmin(tmp,minv[j])) k=j; if(k==-1) break; for(vis[k]=1,sum+=tmp,j=2;j<=n;++j) if(!vis[j]) cmin(minv[j],g[k][j]-f[k][j]*mid); } return sum>=0; } inline double dis(int i,int j){ static double dx,dy; dx=x[i]-x[j],dy=y[i]-y[j]; return sqrt(dx*dx+dy*dy); } int main(){ while(1){ n=read(); if(!n) return 0; for(int i=1;i<=n;++i){ x[i]=read(), y[i]=read(), z[i]=read(); } for(int i=1;i<=n;++i) for(int j=i+1;j<=n;++j) f[j][i]=f[i][j]=dis(i,j), g[j][i]=g[i][j]=abs(z[i]-z[j]); l=0,r=100; while(r-l>eps){ mid=(l+r)/2, (check()?l:r)=mid; } printf("%.3lf\n",l); } return 0; } ``` ### 3.最优比率环 简而言之,最优比率环就是给你一个图,图上每条边有两个权值(当然,第二个权值可能恒为 1 ),然后让你在图中找到一个环,令**环上边的个两个权值和的比值**最大(或是最小) 然后就是要用上面生成树类似的思路,得到每条边边权为: $a[i]-b[i]\times ans$ ,然后在图上找负环(有时是正环,但是正负环可以人工转化的嘛),找到就 $check$ 成功 #### 题目 [[HNOI2009]最小圈](https://www.luogu.org/problemnew/show/P3199) #### 代码 ``` //by Judge #include<cstdio> #include<cstring> #include<iostream> #define eps 1e-10 using namespace std; const int M=1e5+7; inline bool cmin(double& a,double b){return a>b?a=b,1:0;} inline int read(){ int x=0,f=1; char c=getchar(); for(;!isdigit(c);c=getchar()) if(c=='-') f=-1; for(;isdigit(c);c=getchar()) x=x*10+c-'0'; return x*f; } int n,m,tim,pat,head[M],vis[M]; double l,r,mid,d[M]; struct Edge{ int to,next; double val; Edge(){} Edge(int a,double b,int c){ to=a,val=b,next=c; } }e[M<<1]; inline void add(int u,int v,double z){ e[++pat]=Edge(v,z,head[u]),head[u]=pat; } bool dfs(int u){ vis[u]=1; for(int i=head[u];i;i=e[i].next) if(cmin(d[e[i].to],d[u]+e[i].val-mid)) if(vis[e[i].to]||dfs(e[i].to)) return vis[u]=0,1; return vis[u]=0; } inline bool check(){ memset(d,0,sizeof(d)); for(int i=1;i<=n;++i) if(dfs(i)) return 1; return 0; } int main(){ n=read(),m=read(); double z; for(int i=1,x,y;i<=m;++i){ x=read(),y=read(); scanf("%lf",&z); add(x,y,z); if(z<0) l+=z; else r+=z; } while(r-l>eps){ mid=(l+r)/2, (check()?r:l)=mid; } return !printf("%.8lf\n",l); } ``` ### 4.最大密度子图 题目一般就是给出 $n$ 个点然后让你导出一个子图,让这个子图中**边数与点数的比值**最大 不过有时候题目中的边(甚至是点)可能会带上权值,那么具体问题具体分析就好了 #### 题目 [Hard Life](https://www.luogu.org/problemnew/show/UVA1389) (话说这道题翻译照搬的$bzoj$,$UVA$ 是要输出方案的啊...) #### 代码 (好像就这个代码最长了,题目讲得简单然鹅码量巨大) ``` //by Judge #include<cmath> #include<queue> #include<cstdio> #include<cstring> #include<iostream> #define db double using namespace std; const db eps=1e-8;const int N=1005; inline int read(){ int x=0,f=1; char c=getchar(); for(;!isdigit(c);c=getchar()) if(c=='-') f=-1; for(;isdigit(c);c=getchar()) x=x*10+c-'0'; return x*f; } int n,m,S,T,pat,ans,h[N],q[N],head[N],cur[N]; int u[N],v[N],du[N],vis[N]; db U,l,r,mid; struct Edge { int to,nxt; db flow; } e[N<<5]; void insert(int u,int v,db w) { e[++pat]=(Edge){v,head[u],w},head[u]=pat; e[++pat]=(Edge){u,head[v],0},head[v]=pat; } #define v e[i].to inline bool bfs() { int hd=0,tl=1,u; memset(h,0,sizeof(h)),h[q[1]=S]=1; while(hd<tl) { u=q[++hd]; for(int i=head[u]; i; i=e[i].nxt) if(e[i].flow>eps&&!h[v]) h[v]=h[u]+1,q[++tl]=v; } return h[T]; } db dfs(int x,db F) { if(x==T) return F; db w,used=0; for(int &i=cur[x]; i; i=e[i].nxt) if(h[x]+1==h[v]) { w=dfs(v,min(e[i].flow,F-used)); e[i].flow-=w,e[i^1].flow+=w; used+=w; if(used==F) return F; } if(!used) h[x]=0; return used; } db dinic() { db res=0; for(;bfs();res+=dfs(S,1e16)) memcpy(cur,head,sizeof(cur)); return res; } void DFS(int x) { ++ans,vis[x]=1; for(int i=head[x]; i; i=e[i].nxt) if(e[i].flow>eps&&!vis[v]) DFS(v); } #undef v inline bool check() { pat=1,S=0,T=n+1; memset(head,0,sizeof(head)); for(int i=1; i<=n; ++i){ insert(S,i,U), insert(i,T,U+mid+mid-du[i]); } for(int i=1; i<=m; ++i){ insert(u[i],v[i],1), insert(v[i],u[i],1); } return U*n-dinic()>eps; } int main() { while(scanf("%d%d",&n,&m)!=EOF) { for(int i=0; i<=n; ++i) du[i]=0; memset(vis,0,sizeof(vis)),ans=-1,U=m; if(!m) {puts("1\n1\n\n");continue;} for(int i=1; i<=m; ++i) u[i]=read(),v[i]=read(), ++du[u[i]],++du[v[i]]; for(l=0,r=m;r-l>1.0/n/n;){ mid=(l+r)/2, (check()?l:r)=mid; } mid=l,check(),DFS(S); printf("%d\n",ans); for(int i=1; i<=n; ++i) if(vis[i]) printf("%d\n",i); } return 0; } ``` ### 5.最优比率流 还是那句老话,和之前一样的,就是根据二分出来的答案赋边权,然后一般是每次建一边图跑费用流根据最小费用的正负性判断 $ans$ 合法性(用 $zkw$ 跑费用流有时候会挂,比如这个例题,当然可能是我天生大常数那就 $GG$ ) #### 题目 [[SDOI2017]新生舞会](https://www.luogu.org/problemnew/show/P3705) (这里是一道二分图最大费用最大流,但是代码实现中可以边权取反然后跑最小费用) #### 代码 ``` //by Judge #include<queue> #include<cstdio> #include<cstring> #include<iostream> #define eps 1e-8 #define ll long long using namespace std; const double inf=1e18+7; const int M=205; #ifndef Judge #define getchar() (p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<21,stdin),p1==p2)?EOF:*p1++) #endif inline int Min(int a,int b){return a<b?a:b;} inline bool cmax(double& a,double b){return a<b?a=b,1:0;} char buf[1<<21],*p1=buf,*p2=buf; inline int read(){ int x=0,f=1; char c=getchar(); for(;!isdigit(c);c=getchar()) if(c=='-') f=-1; for(;isdigit(c);c=getchar()) x=x*10+c-'0'; return x*f; } int n,S,T,pat=1,tot,head[M],cur[M]; ll a[M][M],b[M][M]; int inq[M],vis[M]; double l,r=1e4,mid,ans,dis[M]; queue<int> q; struct Edge{ int to,flow,next; double val; Edge(int a,double b,int c,int d){ to=a,val=b,flow=c,next=d; } Edge(){} }e[M*M*M]; inline void add(int u,int v,double c,int w){ e[++pat]=Edge(v,c,w,head[u]),head[u]=pat; e[++pat]=Edge(u,-c,0,head[v]),head[v]=pat; } #define v e[i].to inline bool spfa(){ for(int i=S;i<=T;++i) dis[i]=-inf,vis[i]=0; dis[S]=0,q.push(S); while(!q.empty()){ int u=q.front(); q.pop(),vis[u]=0; for(int i=head[u];i;i=e[i].next) if(e[i].flow&&cmax(dis[v],dis[u]+e[i].val)) if(!vis[v]) vis[v]=1,q.push(v); } return dis[T]!=-inf; } int dfs(int u,int flow){ if(u==T) return ans+=dis[T]*flow,flow; int res=0; vis[u]=1; for(int &i=cur[u],k;i;i=e[i].next) if(dis[v]==dis[u]+e[i].val) if(!vis[v]&&e[i].flow){ k=dfs(v,Min(e[i].flow,flow-res)); if(k){ res+=k; e[i].flow-=k; e[i^1].flow+=k; if(res==flow) break; } } return res; } inline bool check(){ pat=1,ans=0, memset(head,0,sizeof(head)); for(int i=1;i<=n;++i){ add(S,i,0,1), add(n+i,T,0,1); } for(int i=1;i<=n;++i) for(int j=1;j<=n;++j) add(i,n+j,a[i][j]-b[i][j]*mid,1); for(;spfa();dfs(S,1e9)) memcpy(cur,head,sizeof(cur)); return ans<=0; } int main(){ n=read(),T=n+1<<1; for(int i=1;i<=n;++i) for(int j=1;j<=n;++j) a[i][j]=read(); for(int i=1;i<=n;++i) for(int j=1;j<=n;++j) b[i][j]=read(); for(tot=pat;r-l>eps;){ mid=(l+r)/2, (check()?r:l)=mid; } return !printf("%.6lf\n",l); } ``` ### 6.各种其他结合 ## ① 与树上背包结合 题意:给你一棵 $n$ 个节点的树,每个节点都有体积和价值。要求选出一些点,使得这些点的体积不小于 $W$ ,且总价值除以总体积最大,求这个最大值。 讲道理就是最优比率背包转移到了树上,没什么特别的(~~就是特别的烦~~)。 题目:咱出到[牛客](https://ac.nowcoder.com/acm/contest/2271/F)上了 题解:也在[牛客](https://ac.nowcoder.com/discuss/346481?type=101&order=0&pos=3&page=1)上 #### 变形: 1. 存在限制条件:相邻的点不可同时选 2. 存在限制条件:一个节点被选择时它的父节点必须被选择(即依赖性背包) 3. 存在限制条件:与一个被选择节点直接相连被选择节点不得超过 k 个 ## 分数规划的一点感想 讲道理,其实分数规划这个东西哪儿都能套,然鹅现在关于分数规划的题目却并不是很多,至少没有到泛滥的地步,不过相信在不(知道多)久的将来,分数规划会成为人尽皆知的一种算法 $$tHe\ EnD$$ 哦,对了,参考文献? [01分数规划问题相关算法与题目讲解](https://blog.csdn.net/hzoi_ztx/article/details/54898323) [poj 3155 (最大密度子图)](https://www.xuebuyuan.com/3225634.html) 另外的另外,**有时间的话**会把自己出的分数规划的题目放到**6**这块内容上来