[Algorithm]01分數規划——Update:2012年7月27日


【關鍵字】

0/1分數規划、最優比率生成樹、最優比率環

【背景】

 根據樓教主的回憶錄,他曾經在某一場比賽中秒掉了一道最優比率生成樹問題,導致很多人跟風失敗,最終悲劇。

自己總結了一些這種問題的解法,因為水平有限,如果有錯誤或是麻煩的地方,盡管噴,郵箱或是下方留言。

聯系我的話perseawe@163.com,歡迎討論,請在標題前注明[acm]或是[oi],以免被垃圾郵件。

 

【知識儲備】

只會用到簡單的公式的整理與變形,還有求和sigma。

 

【定義】

01分數規划問題:所謂的01分數規划問題就是指這樣的一類問題,給定兩個數組,a[i]表示選取i的收益,b[i]表示選取i的代價。如果選取i,定義x[i]=1否則x[i]=0。每一個物品只有選或者不選兩種方案,求一個選擇方案使得R=sigma(a[i]*x[i])/sigma(b[i]*x[i])取得最值,即所有選擇物品的總收益/總代價的值最大或是最小。

01分數規划問題主要包含一般的01分數規划、最優比率生成樹問題、最優比率環問題、最大密度子圖等。我們將會對這四個問題進行討論。

 

永遠要記得,我們的目標是使R取到最值。這句話我會在文中反復的強調。

 

【一些分析】

 數學分析中一個很重要的方法就是分析目標式,這樣我們來看目標式。

R=sigma(a[i]*x[i])/sigma(b[i]*x[i])

我們來分析一下他有什么性質可以給我們使用。

我們先定義一個函數F(L):=sigma(a[i]*x[i])-L*sigma(b[i]*x[i]),顯然這只是對目標式的一個簡單的變形。分離參數,得到F(L):=sigma((a[i]-L*b[i])*x[i])。這時我們就會發現,如果L已知的話,a[i]-L*b[i]就是已知的,當然x[i]是未知的。記d[i]=a[i]-L*b[i],那么F(L):=sigma(d[i]*x[i]),多么簡潔的式子。我們就對這些東西下手了。

再次提醒一下,我們的目標是使R取到最大值。

我們來分析一下這個函數,它與目標式的關系非常的密切,L就是目標式中的R,最大化R也就是最大化L。

F的值是由兩個變量共同決定的,即方案X和參數L。對於一個確定的參數L來說,方案的不同會導致對應的F值的不同,那么這些東西對我們有什么用呢?

假設我們已知在存在一個方案X使得F(L)>0,這能夠證明什么?

F(L)=sigma(a[i]*x[i])-L*sigma(b[i]*x[i])>0即sigma(a[i]*x[i])/sigma(b[i]*x[i])>L也就是說,如果一個方案使得F(L)>0說明了這組方案可以得到一個比現在的L更優的一個L,既然有一個更優的解,那么為什么不用呢?

顯然,d數組是隨着L的增大而單調減的。也就是說,存在一個臨界的L使得不存在一種方案,能夠使F(L)>0. 我們猜想,這個時候的L就是我們要求的最優解。之后更大的L值則會造成無論任何一種方案,都會使F(L)<0.類似於上面的那個變形,我們知道,F(L)<0是沒有意義的,因為這時候的L是不能夠被取得的。當F(L)=0使,對應方案的R值恰好等於此時的L值。

 綜上,函數F(L)有這樣的一個性質:在前一段L中可以找到一組對應的X使得F(L)>0,這就提供了一種證據,即有一個比現在的L更優的解,而在某個L值使,存在一組解使得F(L)=0,且其他的F(L)<0,這時的L無法繼續增大,即這個L就是我們期望的最優解,之后的L會使得無論哪種方案都會造成F(L)<0.而我們已經知道,F(L)<0是沒有任何意義的,因為此時的L值根本取不到。

最后一次提醒,我們的目標是R!!!

如果現在你覺得有些暈的話,那么我要提醒你的就是,千萬不要把F值同R值混淆。F值是根據我們的變形式求的D數組來計算的,而R值則是我們所需要的真實值,他的計算是有目標式決定的。F值只是提供了一個證據,告訴我們真正最優的R值在哪里,他與R值本身並沒有什么必然的聯系。

根據這樣的一段性質,很自然的就可以想到二分L值,然后驗證是否存在一組解使得F(L)>0,有就移動下界,沒有就移動上界。

 所有的01分數規划都可以這么做,唯一的區別就在於求解時的不同——因為每一道題的限制條件不同,並不是每一個解都是可行解的。比如在普通的數組中,你可以選取1、2、3號元素,但在生成樹問題中,假設1、2、3號元素恰好構成了一個環,那就不能夠同時選擇了,這就是需要具體問題,具體分析的部分。

二分是一個非常通用的辦法,但是我們來考慮這樣的一個問題,二分的時候我們只是用到了F(L)>0這個條件,而對於使得F(L)>0的這組解所求到的R值沒有使用。因為F(L)>0,我們已經知道了R是一個更優的解,與其漫無目的的二分,為什么不將解移動到R上去呢?求01分數規划的另一個方法就是Dinkelbach算法,他就是基於這樣的一個思想,他並不會去二分答案,而是先隨便給定一個答案,然后根據更優的解不斷移動答案,逼近最優解。由於他對每次判定使用的更加充分,所以它比二分會快上很多。但是,他的弊端就是需要保存這個解,而我們知道,有時候驗證一個解和求得一個解的復雜度是不同的。二分和Dinkelbach算法寫法都非常簡單,各有長處,大家要根據題目謹慎使用。

 

【實踐】

上面啰嗦了這么多,現在給出程序的框架。

二分法

L:=...;R:=...;
Repeat
  Mid:=(L+R)/ 2;
  For I= 1..X  do D[i]:=A[i]-Mid*B[i];//根據Mid計算D數組
   if 檢查(Mid)成功  then L:=Mid  else R:=Mid;

Until abs(L-R)<Eps;   

Dinkelbach算法
L:=隨便什么東西;
Repeat
  Ans:=L;
  For I= 1..X  do D[i]:=A[i]-L*B[i];//根據L計算D數組
  檢查解並記錄;
  p:= 0;q:= 0;
   for I=每一個元素  do 
     如果元素I在解中
         begin
          p:=p+A[i];q:=q+B[i];
         end;
  L:=p/q;//更新解
Until abs(Ans-L)<Eps;

 

其中檢查解的部分是要看具體情況的。

 

【例題1Poj2976Dropping tests——普通的01分數規划】

大意:給定A數組B數組,從中選擇N-K個使得R最大,輸出Round(100*R);

分析:限制很簡單,只是數目上有所限制,處理方法也很簡單,求出D數組后從大到小排列,從先前向后取N-K個即可,這時的D一定是最大的。

代碼:

二分代碼 110MS
 1 // 10147353      perseawe         2976    Accepted        896K    110MS   Pascal  1517B    2012- 05- 03  10: 09: 47
 2 
 3 Const
 4   Eps= 1e-6;
 5 
 6 Var
 7   n,k:Longint;
 8   Ans:Double;
 9   a,b,c:Array [ 0.. 1000+ 100of Longint;
10   d:Array [ 0.. 1000+ 100of Double;
11 
12 Procedure Init;
13    var
14     i:longint;
15    begin
16     readln(n,k);
17      if (n= 0) and(k= 0then Halt;
18      for i:= 1  to n  do read(a[i]);readln;
19      for i:= 1  to n  do read(b[i]);readln;
20    end;
21 
22  procedure swap( var a,b:Longint); var t:Longint; begin t:=a;a:=b;b:=t; end;
23  procedure swap( var a,b:double); var t:Double; begin t:=a;a:=b;b:=t; end;
24 
25 Procedure Qsort(l,r:Longint);
26    var
27     a,b:Longint;
28     mid:Double;
29    begin
30     a:=l;b:=r;mid:=d[(l+r) shr  1];
31      repeat
32        while d[a]>mid  do inc(a);
33        while d[b]<mid  do dec(b);
34        if a<=b  then
35          begin
36           swap(d[a],d[b]);
37           swap(c[a],c[b]);
38           inc(a);dec(b);
39          end;
40      until a>=b;
41      if a<r  then qsort(a,r);
42      if l<b  then qsort(l,b);
43    end;
44 
45 Procedure Main;
46    var
47     m,i:Longint;
48     L,R,Mid,tmp:Double;
49    begin
50     // 2 '
51      m:=n-k;
52 
53     Mid:= 0;
54      for i:= 1  to n  do  if a[i]/b[i]>Mid  then Mid:=a[i]/b[i];
55     L:= 0;R:=Mid;
56 
57     Repeat
58       Mid:=(L+R)/ 2;
59        for i:= 1  to n  do  begin d[i]:=a[i]-Mid*b[i];c[i]:=i; end;
60       Qsort( 1,n);
61       tmp:= 0;
62        for i:= 1  to m  do tmp:=tmp+d[i];
63        if tmp> 0  then L:=Mid  else R:=Mid;
64     Until abs(L-R)<Eps;
65     Ans:=L;
66    end;
67 
68 Procedure Print;
69    begin
70     writeln(Round(Ans* 100));
71    end;
72 
73 Begin
74   While True Do
75      begin
76       Init;
77       Main;
78       Print;
79      end;
80 End.
Dinkelbach代碼 32MS
 1 // 10147329      perseawe         2976    Accepted        896K    32MS    Pascal  1455B    2012- 05- 03  10: 02: 32
 2 
 3 Const
 4   Eps= 1e-6;
 5 
 6 Var
 7   n,k:Longint;
 8   Ans:Double;
 9   a,b,c:Array [ 0.. 1000+ 100of Longint;
10   d:Array [ 0.. 1000+ 100of Double;
11 
12 Procedure Init;
13    var
14     i:longint;
15    begin
16     readln(n,k);
17      if (n= 0) and(k= 0then Halt;
18      for i:= 1  to n  do read(a[i]);readln;
19      for i:= 1  to n  do read(b[i]);readln;
20    end;
21 
22  procedure swap( var a,b:Longint); var t:Longint; begin t:=a;a:=b;b:=t; end;
23  procedure swap( var a,b:double); var t:Double; begin t:=a;a:=b;b:=t; end;
24 
25 Procedure Qsort(l,r:Longint);
26    var
27     a,b:Longint;
28     mid:Double;
29    begin
30     a:=l;b:=r;mid:=d[(l+r) shr  1];
31      repeat
32        while d[a]>mid  do inc(a);
33        while d[b]<mid  do dec(b);
34        if a<=b  then
35          begin
36           swap(d[a],d[b]);
37           swap(c[a],c[b]);
38           inc(a);dec(b);
39          end;
40      until a>=b;
41      if a<r  then qsort(a,r);
42      if l<b  then qsort(l,b);
43    end;
44 
45 Procedure Main;
46    var
47     p,q:Int64;
48     m,i:Longint;
49     L:Double;
50    begin
51     //Dinkelbach
52     m:=n-k;
53     l:= 1;
54     Repeat
55       Ans:=L;
56        for i:= 1  to n  do  begin d[i]:=a[i]-L*b[i];c[i]:=i; end;
57       Qsort( 1,n);
58       p:= 0;q:= 0;
59        for i:= 1  to m  do
60          begin
61           inc(p,a[c[i]]);
62           inc(q,b[c[i]]);
63          end;
64       L:=p/q;
65     Until abs(L-Ans)<Eps;
66    end;
67 
68 Procedure Print;
69    begin
70     writeln(Round(Ans* 100));
71    end;
72 
73 Begin
74   While True Do
75      begin
76       Init;
77       Main;
78       Print;
79      end;
80 End.

另外:如果是最小選擇N-K個怎么辦?       

辦法是一樣的,從大到小排列序,傻子才多選,能少選就少選。反正F值具體的大小沒什么關系,我們只要知道他與0的關系即可。

 

【例題2Poj2728Desert King——最優比率生成樹】

大意:給定一張圖,每條邊有一個收益值和一個花費值,求一個生成樹,要求花費/收益最小,輸出這個值

分析:現在的限制就有點復雜了,要求解必須是一棵生成樹。而且這道題目要求的花費/收益最小,當然你求收益/花費最大然后反過來也是可以的,注意處理花費為0的情況。如果求最小的,處理方法是也類似的,先求個D,然后做一次最小生成樹,顯然得到的就是函數值。不過這道題用Dinkelbach比二分好的多。

Dinkelbach代碼
 1 // 10148420      perseawe         2728    Accepted        916K    407MS   Pascal  1560B    2012- 05- 03  16: 03: 10
 2 
 3 Const
 4   Eps= 1e-6;
 5   MaxN= 1000+ 100;
 6 
 7 Var
 8   n:Longint;
 9   ans:Double;
10   x,y,h:Array [ 0..MaxN]  of Longint;
11   Use:Array [ 0..MaxN]  of Boolean;
12   a,b,d:Array [ 0..MaxN]  of Double;
13 
14 Procedure Init;
15    var
16     i:Longint;
17    begin
18     readln(n);
19      if n= 0  then Halt;
20      for i:= 1  to n  do readln(x[i],y[i],h[i]);
21    end;
22 
23 Procedure Main;
24    var
25     i,m,pos:Longint;
26     L,tmp,ta,tb,p,q:Double;
27    begin
28     L:= 0;
29     Repeat
30       Ans:=L;
31       //Prim
32       Fillchar(Use,sizeof(Use),False);Use[ 1]:=True;
33       For i:= 2  to n  do
34          begin
35           a[i]:=abs(h[i]-h[ 1]);
36           b[i]:=sqrt(sqr(x[i]-x[ 1])+sqr(y[i]-y[ 1]));
37           d[i]:=a[i]-L*b[i];
38          end;
39       m:= 1;p:= 0;q:= 0;
40       While m<n  do
41          begin
42           tmp:= 1000000000;
43            for i:= 2  to n  do
44              if  not(Use[i]) and(d[i]<tmp)  then
45                begin
46                 tmp:=d[i];pos:=i;
47                end;
48           Use[pos]:=true;p:=p+a[pos];q:=q+b[pos];
49            for i:= 2  to n  do
50              if  not(Use[i])  then
51                begin
52                 ta:=abs(h[i]-h[pos]);tb:=sqrt(sqr(x[i]-x[pos])+sqr(y[i]-y[pos]));
53                  if ta-L*tb<d[i]  then
54                    begin
55                     d[i]:=ta-L*tb;
56                     a[i]:=ta;b[i]:=tb;
57                    end;
58                end;
59           Inc(m);
60          end;
61       L:=p/q;
62     Until abs(L-Ans)<Eps;
63    end;
64 
65 Procedure Print;
66    begin
67     writeln(ans: 0: 3);
68    end;
69 
70 Begin
71    while True Do
72      begin
73       Init;
74       Main;
75       Print;
76      end;
77 End.

 最小生成樹用了Prim,只要不是實在沒辦法,還是不要在稠密圖特別是完全圖上用Kruskal。

 

【例題3Poj3621Sightseeing Cows——最優比率環】

大意:給定一張圖,邊上有花費,點上有收益,點可以多次經過,但是收益不疊加,邊也可以多次經過,但是費用疊加。求一個環使得收益和/花費和最大,輸出這個比值。

分析:比上面更加的惡心了。先不說環的問題,就是花費和收益不在一處也令人蛋疼。這時候需要用到幾個轉化和結論。

     首先的一個結論就是,不會存在環套環的問題,即最優的方案一定是一個單獨的環,而不是大環套着小環的形式。這個的證明其實非常的簡單,大家可以自己想一下(提示,將大環上的收益和記為x1,花費為y1,小環上的為x2,y2。重疊部分的花費為S。表示出來分類討論即可)。有了這個結論,我們就可以將花費和收益都轉移到邊上來了,因為答案最終一定是一個環,所以我們將每一條邊的收益規定為其終點的收益,這樣一個環上所有的花費和收益都能夠被正確的統計。

     解決了蛋疼的問題之后,就是01分數規划的部分了,我們只需要計算出D數組后找找有沒有正權環即可,不過這樣不太好,不是我們熟悉的問題,將D數組全部取反之后,問題轉換為查找有沒有負權環,用spfa或是bellman_ford都可以。這道題目就是典型的不適合用Dinkelbach,記錄一個負權環還是比較麻煩的,所以二分搞定。

二分代碼
 1 // 10148804      perseawe         3621    Accepted        1000K   422MS   Pascal  1239B    2012- 05- 03  17: 17: 40
 2 
 3 Const
 4   Eps= 1e-6;
 5   MaxNode= 1000+ 100;
 6   MaxEdge= 5000+ 500;
 7 
 8 Var
 9   Ans:Double;
10   n,m:Longint;
11   a:Array [ 0..MaxNode]  of Longint;
12   dis: array [ 0..MaxNode]  of Double;
13   b,u,v:Array [ 0..MaxEdge]  of Longint;
14   d:Array [ 0..MaxEdge]  of Double;
15 
16 Procedure Init;
17    var
18     i:Longint;
19    begin
20     readln(n,m);
21      for i:= 1  to n  do readln(a[i]);
22      for i:= 1  to m  do readln(u[i],v[i],b[i]);
23    end;
24 
25 Function Bellman_Ford(L:Double):Boolean;
26    var
27     i,j:Longint;
28     Flag:Boolean;
29    begin
30      for i:= 1  to m  do d[i]:=-(a[v[i]]-L*b[i]);
31      for i:= 1  to n  do dis[i]:= 0;
32      for i:= 1  to n  do
33        begin
34         Flag:=False;
35          for j:= 1  to m  do
36            if Dis[u[j]]+d[j]<Dis[v[j]]  then
37              begin
38               Dis[v[j]]:=Dis[u[j]]+d[j];
39               Flag:=True;
40              end;
41         If  not(Flag)  then Exit(False);
42        end;
43     Exit(True);
44    end;
45 
46 Procedure Main;
47    var
48     L,R,Mid:Double;
49    begin
50     L:= 0;R:= 20000;
51     Repeat
52       Mid:=(L+R)/ 2;
53        if Bellman_Ford(Mid)  then L:=Mid  else R:=Mid;
54     Until abs(L-R)<Eps;
55     Ans:=L;
56    end;
57 
58 Procedure Print;
59    begin
60      if Ans>Eps  then writeln(Ans: 0: 2else writeln( 0);
61    end;
62 
63 Begin
64   Init;
65   Main;
66   Print;
67 End.

因為圖省事,所以用了Bellman_Ford。還有就是注意無解的判斷,無解時檢查會一直不成功的,所以上界會一直向下移動的。

 

【后記】

本來Zoj上還有一道題的,時間不夠了就先放過去吧, 但是常見的三種01規划我都已經列舉到了並給出了常見的處理手段。

算法運用之妙,存乎一心也。數學是最神奇的。強烈建議大家學好數學!

本來還有一個例題0和一個非常巧妙的數學證明的,但是發覺沒有什么太大的意義,就省略掉了。大意就是給定A數組和B數組(A、B的元素值都大於0),最小取一個,求最大的R值。顯然直接計算所有的A[i]/B[i]取最大值就可以,因為兩個分數分子分母對應相加得到的結果一定是小於較大的那個。具體的證明很簡單,分類討論即可。

另一個被我省略的部分就是對Dinkelbach算法的分析,這需要更強的數學分析才行,因為這並不是重點,所以我將其省略掉了。如果有興趣的同學可以找一下一篇叫做《對於0-1分數規划的Dinkelbach算法的分析》,由武鋼三中 吳豪[譯]的一篇文章看一下。

最后感謝網上很多大牛的題解和心得,特別是This_poet的題解,第一二題中參考了她的代碼,幾乎已經是Copy了,見諒。>_<

歡迎大家拍磚討論,具體看背景。

 

————————————————Update:2012年5月10日——————————————

【例題4Poj3155——Hard Life最大密度子圖】

大意:給定G=(V,E)求其中的一個子圖使得邊數/點數最大

分析:詳見《最小割模型在信息學競賽中的應用》作者胡伯濤。

 

 ————————————————Update:2012年6月19日——————————————

【例題5Zoj2676——Network Wars】

大意:給定一張圖,規定一個割的平均值是邊權和/邊數.求平均值最小的割.

分析:非常顯然,因為是求最小,所以只要對於某個L值,解空間中最小的一組解滿足g(x)=0即可。於是,二分L后改變權值求解最小割。但是會出現負權的情況,這里需要特殊處理。負權是一定會出現在解中的,遇到負權直接加上即可。如果最小割<0則L增大,反之減少.還有就是注意精度的問題,小心處理。

 【例題6游戲——最大密度子圖變種】

大意:給定N個人,可以選擇任意多個人,記為K。單個人是沒有戰斗力的,必須要合體才有戰斗力(大霧),給定所有的Aij表示當選擇方案同時出現i和j兩個人時的戰斗力加成。特殊的Aii=0。但是選擇人不是沒有代價的,代價是k*(200-k).定義這個方案的評分=sigma(A[i,j]|if i,j both selected)/(k*(200-k))。要求最大化評分。

 N<=50...也就是滿足選擇的人越多代價越多.

分析:顯然是一個最大密度子圖的模型。但是與一般的模型中,表達式中的分母不僅僅含有點數,更含有一個點數的平方。這是令人非常難受的。問題的核心就在於如何將這個平方蘊含到圖里面去。觀察到這個圖非常的特殊,當選擇一些點時,這些點構成的子圖是一個完全圖,即邊數是K*(K-1)/2.哈哈,這里也有一個K的平方,於是問題得以解決了。二分L后整理表達式,將分母上的K的平方蘊含到邊權上去。這樣就可以轉化為一般的最大密度子圖。

【一些理解】

從這篇文章寫出來到現在也已經很長時間了,對01分數規划問題也有了很多新的看法。

以下是一些個人的理解,很有可能存在錯誤。大家幫忙看看有沒有問題。

01分數規划問題求最值,但是最值有兩種,一種是最大值,一種是最小值。

在我看來,對於每一個L。我們先假定L是確定了的,這時問題還會有很多的方案,每一個方案有其評估值。我們求最小值時是找這樣的一個L:所有的合法方案中,只有唯一的一個方案評估值為0,其他的方案評估值都>0.而最大值恰好相反,即也是只有唯一的一個方案評估值為0,但是其他的方案評估值都<0.

在上面的推導中,我已經說過了,當存在評估值>0的方案時,說明L是可以增大的。因為計算方案對應的原始表達式值計算出的L'一定是大於L的,同理,而評估值<0的方案則是不優的。

雖然題目要求的東西會不同,生成樹呀,割呀,子圖之類的呀。但是我覺得上面的這個東西是一個通用的玩意,是01分數規划本身決定的,而不是題目決定的。我現在還是有一個概念不是特別明白。先寫到這里,歡迎討論。

————————————————Update:2012年7月27日————————————————————————————

 

 【例題7Poj3266——CowSchool0/1分數規划+數據結構】
大意:(From applepi)
Bessie考了N科(N≤50000),每科的得分為Ti,滿分為Pi(樣例似乎是五科紅燈……)。現在Bessie的老師要給她們統計最終得分,方法 是先算出每一科得分的百分比,去除D個百分比最低的科目,然后剩下科目的∑T / ∑P就是最終得分的百分比。很幸運,沒有人有兩科的得分比一樣。Bessie的數學很牛X,她馬上就發現她的老師在坑爹,因為有時候可以通過去除不同的科 目來獲得更高的分數。現在Bessie想知道,對於哪些D值,她可以通過去除與老師不同的D科,從而獲得更高的分數呢?
題解:顯然,Bessie的老師沒有看過這篇文章(>_<)。要不他(她)就不會犯這么意識流的問題了。現在題目要求我們給出一個"證據"來證明那些D老師的做法並不是最優的.由上面的一些結論,我們可以知道,只要我們能夠找到在老師給定的百分比下的一組解,其值>0即可,這就證明了存在一個百分比要優於老師的解.這時我們可以用一個巧妙的轉化,顯然,老師選定的解在加權之后的和一定是0.那么,一個簡單的構造方法是找出老師選定中的最小加權值以及沒選定的中最大的加權值,如果后者大於前者,那么顯然將前者更換為后者得到的解值>0,說明存在更優的解。然后這樣的一個問題可以把每一個分數看做平面上的一個點,用平衡樹維護凸殼或是充分挖掘單調性后使用棧和隊列來維護都是可以快速的完成查找操作的。具體見applepi的題解就好。

 

這道題就把01分數規划作為一個工具,要求挖掘更加符合這個題目的特性,還是很難的。說實話,這道題是Poj上Usaco題目中最難的一道題目了。


免責聲明!

本站轉載的文章為個人學習借鑒使用,本站對版權不負任何法律責任。如果侵犯了您的隱私權益,請聯系本站郵箱yoyou2525@163.com刪除。



 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM