集合上的动态规划

例题：最优配对问题

　　空间里有n个点P₀,P₁,...,P_n-1，你的任务是把他们配成n/2对（n是偶数），使得每个点恰好在一个点对中。所有点对中两点的距离之和应尽量小。n<=20,|x_i|,|y_i|,|z_i|<=1000.

样例输入：

20
1 2 3
1 1 1
5 6 2
4 7 8
2 3 1
1 4 7
2 5 8
3 6 9
1 2 5
2 3 6
4 5 2
7 8 5
4 5 1
-1 2 3
-1 -9 -7
0 0 0
100 0 0
9 5 1
7 5 3
5 5 5

分析：

　　既然每个点都要配对，很容易把问题看成如下多阶段决策过程：先确定P₀和谁配对，然后是P₁，接下来是P₂，……，最后是P_n-1.按照前面的思路，设d(i)表示把前i个点两两配对的最小距离和，然后考虑第i个点的决策——它和谁配对呢？假设它和点j配对（j< i），那么接下来的问题应是“把前i-1个点中除了j之外的其它点两两配对”，它显然无法用任何一个d值来刻画——我们的状态定义无法体现出“除了一些点之外”这样的限制。

　　当发现状态无法转移后，常见的方法是增加维度，即增加新的因素，更细致地描述状态。既然刚才提到了“除了某些元素之外”，不妨把它作为状态的一部分，设d(i,S)表示把前i个点中，位于集合S中的元素两两配对的最小距离，则状态转移方程为：

　　d(i,S) = min{|P_iP_j|+ d(i-1,S-{i}-{j}) | j属于S}

　　其中|PiPj表示Pi和Pj之间的距离。边界是d(-1,S )= 0.

记忆化搜索的代码如下：

 1 #include<cstdio>
 2 #include<cstdlib>
 3 #include<cmath>
 4 #define maxn 20
 5 #define INF (1<<30)
 6 int n;
 7 struct Point{
 8     double x, y, z;
 9 } p[maxn];
10 double f[(1<<maxn)];
11 inline double dist(int i, int j)
12 {
13     return sqrt((p[i].x-p[j].x)*(p[i].x-p[j].x) + (p[i].y-p[j].y)*(p[i].y-p[j].y) + (p[i].z-p[j].z)*(p[i].z-p[j].z));
14 }
15 inline double min(double x, double y)
16 {
17     if(x < y) return x; 
18     return y;
19 }
20 double d(int S){
21     //记忆化搜索快得要逆天了！！（因为奇数的都没有算！！ ）
22     if(f[S] != -1) return f[S];
23     int i, j;
24     double& ans = f[S];
25     ans = (S == 0? 0 : INF);
26     for(i = n-1; i >= 0; --i) if(S & (1<<i)) break;
27     for(j = 0; j < i; ++j) if(S & (1<<j))
28       ans = min(ans, d(S^(1<<i)^(1<<j)) + dist(i, j));
29     return ans;
30 }
31 
32 int main(){
33     int i;
34     scanf("%d", &n);
35     if(n%2) {printf("No answer!
"); return 0;}
36     for(i = 0; i < n; ++i) scanf("%lf%lf%lf", &p[i].x, &p[i].y, &p[i].z);
37     for(i = 0; i < (1 << n); ++i) f[i] = -1;
38     printf("%.3lf
", d((1<<n)-1));
39     /*
40     int cnt = 0;
41     for(i = 0; i < (1 << n); ++i) if(f[i]!=-1) ++cnt;
42     printf("%d
", cnt);
43     */
44     return 0;
45 }

动态规划的代码如下：

 1 #include<cstdio>
 2 #include<cstdlib>
 3 #include<cmath>
 4 #define maxn 20
 5 #define INF (1<<30)
 6 int n;
 7 struct Point{
 8     double x, y, z; 
 9 } p[maxn];
10 double f[maxn][(1<<maxn)];   
11 inline double dist(int i, int j)
12 {
13     return sqrt((p[i].x-p[j].x)*(p[i].x-p[j].x) + (p[i].y-p[j].y)*(p[i].y-p[j].y) + (p[i].z-p[j].z)*(p[i].z-p[j].z));
14 }
15 inline double min(double x, double y)
16 {
17     if(x < y) return x; 
18     return y;
19 }
20 int main(){
21     int i,j;
22     scanf("%d", &n);
23     if(n%2) {printf("No answer!
"); return 0;}
24     for(i = 0; i < n; ++i) scanf("%lf%lf%lf", &p[i].x, &p[i].y, &p[i].z);   
25     for(i = 0; i < n; ++i)
26         for(int S = 0; S < (1 << (i+1)); ++S){
27             double& ans = f[i][S];
28             ans = (S==0 ? 0 : INF);
29             if(S & (1 << i)) {
30                 for(j = 0; j < i; ++j) if(S & (1 << j))
31                     ans = min(ans, f[i-1][S^(1<<i)^(1<<j)] + dist(i, j));
32             }
33             else if(i!=0) f[i][S] = f[i-1][S];  
34         }    
35         printf("%.3lf
", f[n-1][(1<<n)-1]);
36         return 0;
37 }

隐含的阶段

　　上面的方程可以进一步简化。事实上，阶段 i 根本不用保存，它已经隐含在S中了——S中最大的元素就是i。这样，可直接用d(S)表示“把S中的元素两两配对的最小距离和”，则状态转移方程为：

　　　　d(S)=min{|P_iP_j|+d(S-{i}-{j}) | j属于S，i=max(S)}

状态有2ⁿ个，每个状态有O(n)种转移方式，总时间复杂度为O(n2ⁿ).

　　值得一提的是，不少人一直在用这样的状态转移方程：

　　　　d(S)=min{|P_iP_j|+d(S-{i}-{j}) | i,j属于S}

　　它和刚才的方程很类似，唯一不同是：i和j都是需要枚举的。这样做虽然也没错，但每个状态的转移次数高达O(n²)，总时间复杂度为O(n²2ⁿ)，比刚才的方法慢。这个例子说明：即使用相同的状态描述，减少决策也很重要。

　　那么如何求出S中的最大元素呢？用一个循环判断即可。当S取遍{0,1,2,...,n-1}的所有子集时，平均判断次数仅为2。

1 for(int S=0; S<(1<<n); S++)
2 {
3     int i,j;
4     d[S] = INF;
5     for(i=0; i<n; i++)
6         if(S && (1<<i)) break;
7     for(j=i+1; j<n; j++)
8         if(S & (1<<j)) d[S] <?= d[S^(1<<i)^(1<<j)] + dist(i, j));
9 }

　　注意，在上述的程序中求出的i是S中的最小元素，而不是最大元素，但这并不影响答案。另外，j的枚举只需从 i+1开始——既然 i 是S中最小元素，则说明其他元素自然均比i大。最后需要说明的是S的枚举顺序。不难发现：如果S‘是S的真子集，则一定有S'<S,因此S递增的顺序计算，需要用到某个d值时，它一定已经算出来了。

　　如果S'是S的真子集，则一定有S'<S。在用递推法实现子集动态规划时，该规则往往可以确定计算部分。

优化后代码如下：

 1 #include<cstdio>
 2 #include<cstdlib>
 3 #include<cmath>
 4 #define maxn 22
 5 #define maxs (1<<maxn)
 6 #define INF 1e10
 7 int n; 
 8 double f[maxs], x[maxn], y[maxn], z[maxn];
 9 inline double min(double x, double y)
10 {
11     return x<y?x:y;
12 }
13 inline double dist(int i, int j)
14 {
15     return sqrt((x[i]-x[j])*(x[i]-x[j]) + (y[i]-y[j])*(y[i]-y[j]) + (z[i]-z[j])*(z[i]-z[j]));
16 }    
17 int main(){
18     scanf("%d", &n);
19     for(int i  = 0; i < n; ++i) scanf("%lf%lf%lf", x+i, y+i, z+i);
20     int U = (1 << n);
21     for(int S = 0; S < U; ++S){
22         int i;
23         for(i = n-1; i >= 0; --i) if(S & (1<<i)) break;
24         f[S] = (S==0) ? 0 : INF;
25         for(int j = 0; j < i ; ++j) if(S & (1<<j)){
26           f[S] = min(f[S], f[S^(1<<i)^(1<<j)] + dist(i, j));
27         }    
28       }
29     printf("%.3lf
", f[U-1]);
30     return 0;
31 }