【POJ】2728 Desert King 最优比率生成树——01分数规划【经典】

最近在刷巨巨们放出来的专题，然后没做几题就卡住了，果然还是太弱了T U T...

这次做到了一题01分数规划求解的生成树问题。

题目大意是这样的：给你一个无向完全图，每条边ｉ都有两个权值，长度a[ i ]，花费b[ i ]，需要选出其中的一些边构造一颗生成树，生成树需要满足条件：∑ b [ i ] / ∑ a [ i ]最小。

这样我还是先来介绍一下01分数规划吧~

给定一个上述的问题，我们可以设R = ∑x[ i ] * b[ i ] / ∑x[ i ] * a[ i ]，其中x[ i ]为0或1表示边i取或者不取。

这里我们先假设对于所有的可行解x[ i ]，a[ i ] * x[ i ]都是正数。

我们需要求R的最小值，那么不妨先设一个子问题Q(L) = ∑x[ i ] * b[ i ] - L * ∑x[ i ] * a[ i ] 的最小值，然后我们需要求满足条件（即它是一棵生成树）的L的最小值。

易知L就是R，所以问题等价于求R的最小值

首先我们先对Q(L)移项得：Q(L) = ∑x[ i ] * ( b[ i ] - R * a[ i ] )。

如果我们已知L，那么b[ i ] - L * a[ i ]也是已知的，所以子问题只和两个参数L，x[ i ]有关。并且对于一个L，x[ i ]的选取不同，将可能会导致Q(L)的取值不同。

现在我们假设Q(L) < 0 ，那么我们能得到什么信息？

根据Q(L) = ∑x[ i ] * ( b[ i ] - L * a[ i ] ) < 0 可以推出∑x[ i ] * b[ i ] / ∑x[ i ] * a[ i ] < L ，那么就说明我们找到了一个比R更优的解！为什么不选择更优一点的解呢~

如果Q(L) < 0，那么这时候比L小的解都是不是最优的，自然不需要了。（其实应该理解为它们都是意义的）

当然如果Q(L) = 0，那么我们找到了一个可行解R=L，但是不是最优解，暂时我们还不知道。

如果函数Q(L)满足单调性，则说明如果存在Q(L) = 0，则有且只有这一个解，并且这个解就是我们需要的答案。

证明：

假设L2 > L1，令X[ i ]是L2的一组解向量，Y[ i ]是L1的一组解向量。

那么Q(L2) = ∑X[ i ] * b[ i ] - L2 * ∑X[ i ] * a[ i ] ， Q(L1) =∑Y[ i ] * b[ i ] - L1 * ∑Y[ i ] * a[ i ]。

易知∑Y[ i ] * b[ i ] - L2 * ∑Y[ i ] * a[ i ] >= ∑X[ i ] * b[ i ] - L2 * ∑X[ i ] * a[ i ](因为X[ i ]，Y[ i ]分别是Q(L2)，Q(L1)最小值对应的解向量，所以Y[ i ]放到L2中去一定不会比最小值还小。

所以有：Q(L2) - Q(L1) = ( ∑X[ i ] * b[ i ] - L2 * ∑X[ i ] * a[ i ] ) -( ∑Y[ i ] * b[ i ] - L1 * ∑Y[ i ] * a[ i ] ) <=

                                           ( ∑Y[ i ] * b[ i ] - L2 * ∑Y[ i ] * a[ i ] ) - ( ∑Y[ i ] * b[ i ] - L1 * ∑Y[ i ] * a[ i ] ) =

                                            ( L1 - L2 ) * ∑Y[ i ] * a[ i ] < 0 ( Y[ i ] > 0，a[ i ] > 0 )

因此我们得到了对于L2 > L1，Q(L2) < Q(L1)，所以函数Q(L)单调递减。

所以根据Q(L)的单调性，我们很容易知道当Q(L) = 0的时候，R = L 取得最小值。

另外，假设L*是Q(L) = 0时的最优解：

Q(L) > 0
当且仅当
L < L*

Q(L) = 0
当且仅当
L = L*

Q(L) < 0
当且仅当
L > L*

现在我们可以通过二分法来寻找R值或者通过Dinkelbach算法迭代求解。

不同的情况下，两种算法会发挥不同的效果，视情况而定。

传送门：【POJ】2728 Desert King

题目大意：给你一个无向完全图，每条边ｉ都有两个权值，长度a[ i ]，花费b[ i ]，需要选出其中的一些边构造一颗生成树，生成树需要满足条件：∑ b [ i ] / ∑ a [ i ]最小。

题目分析：01分数规划——最优比率生成树。经典题！！！

本题二分明显没有迭代的效率高。

二分解法：

求出Q(L) = height[ i ] - R * length[ i ]

Q(L) <= 0 则移动上界，否则移动下界

二分代码：

//1454ms#include <cmath>#include <cstdio>#include <cstring>#include <algorithm>using namespace std ;#define REP( i , n ) for ( int i = 0 ; i < n ; ++ i )#define REPF( i , a , b ) for ( int i = a ; i <= b ; ++ i )#define REPV( i , a , b ) for ( int i = a ; i >= b ; -- i )#define clear( a , x ) memset ( a , x , sizeof a )const int MAXN = 1005 ;const double INF = 1e18 ;const double eps = 1e-5 ;struct Point {double x , y , h ;void input () {scanf ( "%lf%lf%lf" , &x , &y , &h ) ;}} ;struct MST {Point point[MAXN] ;double G[MAXN][MAXN] ;double length[MAXN][MAXN] ;double height[MAXN][MAXN] ;double d[MAXN] ;bool vis[MAXN] ;int fa[MAXN] ;int n ;double dist ( double x , double y ) {return sqrt ( x * x + y * y ) ;}void input () {REP ( i , n )point[i].input () ;REPF ( i , 0 , n - 1 )REPF ( j , i + 1 , n - 1 ) {length[i][j] = length[j][i] = dist ( point[i].x - point[j].x , point[i].y - point[j].y ) ;height[i][j] = height[j][i] = fabs ( point[i].h - point[j].h ) ;}REP ( i , n )height[i][i] = length[i][i] = G[i][i] = 0 ;}void cal ( double r ) {REPF ( i , 0 , n - 1 )REPF ( j , i + 1 , n - 1 )G[i][j] = G[j][i] = height[i][j] - r * length[i][j] ;}double prim ( double r ) {cal ( r ) ;REP ( i , n )d[i] = INF ;clear ( vis , 0 ) ;d[0] = 0 ;fa[0] = 0 ;double ans = 0 , m ;int x ;REP ( _ , n ) {m = INF ;REP ( i , n )if ( !vis[i] && m > d[i] )m = d[x = i] ;vis[x] = 1 ;ans += m ;REP ( i , n )if ( !vis[i] && d[i] > G[x][i] )d[i] = G[x][i] , fa[i] = x ;}return ans ;}void solve () {input () ;double mid , l = 0 , r = 100 ;while ( fabs ( r - l ) > eps ) {mid = ( l + r ) / 2 ;if ( prim ( mid ) <= eps )r = mid ;elsel = mid ;}printf ( "%.3f\n" , l ) ;}} ;MST z ;int main () {while ( ~scanf ( "%d" , &z.n ) && z.n )z.solve () ;return 0 ;}

Dinkelbach算法：

设初值0，第一次带入求出的Q(L)一定大于0，但是得到的L带入第二次得到的一定小于0，然后就不断的迭代直到Q(L) = 0。

迭代代码：

//282ms#include <cmath>#include <cstdio>#include <cstring>#include <algorithm>using namespace std ;#define REP( i , n ) for ( int i = 0 ; i < n ; ++ i )#define REPF( i , a , b ) for ( int i = a ; i <= b ; ++ i )#define REPV( i , a , b ) for ( int i = a ; i >= b ; -- i )#define clear( a , x ) memset ( a , x , sizeof a )const int MAXN = 1005 ;const double INF = 1e18 ;const double eps = 1e-5 ;struct Point {double x , y , h ;void input () {scanf ( "%lf%lf%lf" , &x , &y , &h ) ;}} ;struct MST {Point point[MAXN] ;double G[MAXN][MAXN] ;double length[MAXN][MAXN] ;double height[MAXN][MAXN] ;double d[MAXN] ;bool vis[MAXN] ;int fa[MAXN] ;int n ;double dist ( double x , double y ) {return sqrt ( x * x + y * y ) ;}void input () {REP ( i , n )point[i].input () ;REPF ( i , 0 , n - 1 )REPF ( j , i + 1 , n - 1 ) {length[i][j] = length[j][i] = dist ( point[i].x - point[j].x , point[i].y - point[j].y ) ;height[i][j] = height[j][i] = fabs ( point[i].h - point[j].h ) ;}REP ( i , n )height[i][i] = length[i][i] = G[i][i] = 0 ;}void cal ( double r ) {REPF ( i , 0 , n - 1 )REPF ( j , i + 1 , n - 1 )G[i][j] = G[j][i] = height[i][j] - r * length[i][j] ;}double prim ( double r ) {cal ( r ) ;REP ( i , n )d[i] = INF ;clear ( vis , 0 ) ;d[0] = 0 ;fa[0] = 0 ;double a = 0 , b = 0 , m ;int x ;REP ( _ , n ) {m = INF ;REP ( i , n )if ( !vis[i] && m > d[i] )m = d[x = i] ;a += length[fa[x]][x] ;b += height[fa[x]][x] ;vis[x] = 1 ;REP ( i , n )if ( !vis[i] && d[i] > G[x][i] )d[i] = G[x][i] , fa[i] = x ;}return b / a ;}void solve () {input () ;double res = 0 , tmp ;while ( 1 ) {tmp = prim ( res ) ;if ( fabs ( tmp - res ) <= eps )break ;res = tmp ;}printf ( "%.3f\n" , res ) ;}} ;MST z ;int main () {while ( ~scanf ( "%d" , &z.n ) && z.n )z.solve () ;return 0 ;}