算法导论（九）－－中位数和顺序统计学

在一个由元素组成的集合里，第i个顺序统计量（order statistic）是该集合第i小的元素。例如，最小值是第1个顺序统计量(i=1)，最大值是第n个顺序统计量(i=n)。一个中位数（median）是它所在集合的“中点元素”。当n为奇数时，中位数是唯一的，i=(n+1)/2；当n为偶数时，中位数有两个，一个是i=floor((n+1)/2)（下中位数），另一个是i=ceiling((n+1)/2)（上中位数）。在不考虑n的奇偶性的情况下，中位数总是出现在i=floor((n+1)/2)处和i=ceiling((n+1)/2)处。

最小值和最大值

在一个有n个元素的集合里，要找到最小元素，可以用下面的代码：

MINIMUM(A) {
1      min = A[1];
2         for i = 2 to A.length
3                  if (min > A[i])
4                           min = A[i];
5         return min;
}

为了找到最小元素，每个元素都必须参与比较，所以n-1次比较是必须的。从比较次数来看，算法MINIMUM是最优的。

如果要同时找出最大值和最小值， 我们很容易想到让每一个元素分别和min和max比较，这样每个元素要进行2次比较，而总的比较数就为2*n-2。事实上只需3*floor(n/2)次 比较就可以同时找到最小值和最大值。我们并不让每个元素同时与min和max比较，而是取两个元素先比较，较大的元素与max比较，而较小的元素与min 比较，这样每两个元素要做3次比较，总的比较次数最多为3*floor(n/2)。

以期望线性时间做选择

如 果我们不选最大值或最小值，而是选一个第i小的值。我们可以用一种分治算法——RAMDOMIZED_SELECT，它以快速排序为模型：把数组随机划分 为两部分，A[p...q-1]的元素比A[q]小，A[q+1...r]的元素比A[q]大。与快速排序不同，如果i=q，则A[q]就是要找的第i小 的元素，返回这个值；如果i < q，则说明第i小的元素在A[p...q-1]里；如果i > q，则说明第i小的元素在A[q+1...r]里。

下面是在A[p...r]中找到第i小元素的代码：

RANDOMIZED_SELECT(A, p, r, i) {
1      if p == r
2                  return A[p];
3         q = RANDOMIZED_PARTITION(A, p, r);
4         k = q-p+1;
5         if i == k
6           return A[q];
7         elseif i < k
8                  return RANDOMIZED_SELECT(A, p, q-1, i);
9         else return RANDOMIZED_SELECT(A, q+1, r, i-k);
}

在最坏情况下，数组被划分为n-1和0两部分，而第i个元素总是落在n-1的那部分里，运行时间为Ө(n^2)。但随机算法更关心的是平均情况：

设指示器随机变量X_k = I{子数组A[p...q]中恰有k个元素}
这样，元素数量为n的数组在每次划分时，就会被划分为1...k-1, k, k+1...n三部分。下次递归时，就可以作用在1...k-1（共k-1个元素）区间或k+1...n（共n-k个元素）区间上。考虑X_k所有可能的取值，则可以得到递归式：
T(n) ≤ ∑X_k * T(max(k-1, n-k)) + O(n)
    = ∑(X_k * T(max(k-1, n-k)) + O(n))

由于X_k为1的概率为1/n，取期望值，得到：
E(T(n)) ≤ E[∑(X_k * T(max(k-1, n-k)) + O(n))]
    = ∑E[(X_k * T(max(k-1, n-k))]+ O(n)
    = ∑E[X_k] * E[T(max(k-1, n-k))] + O(n)　（X_k与T(max(k-1, n-k))独立）
    = ∑1/n * E[T(max(k-1, n-k))] + O(n)

考虑max(k-1, n-k)，如果k > ceiling(n/2)，其值为k-1；若k ≤ ceiling(n/2)，其值为n-k。所以
∑T(max(k-1, n-k)) ≤ 2 * ∑T(k)，于是有：
E(T(n)) ≤ 2/n * ∑E[T(k)] + O(n)

用代换法解上式，可以得到E(T(n))=O(n)。也就是说在平均情况下，任何顺序统计量（特别是中位数）都可以在线性时间内得到。

最坏情况线性时间的选择

相比于上面的随机选择，我们有另一种类似的算法，它在最坏情况下也能达到O(n)。它也是基于数组的划分操作，而且利用特殊的手段保证每次划分两边的子数组都比较平衡。

首先我们需要稍微修改一下PARTITION算法（不是RANDOMIZED_PARTITION)，它接收一个数组和一个值x，并把它划分为小于x和大于x的两部分（x为A中某个元素的值）：

PARTITION_X(A, p, r, x) {
1      for i = 1 to n
2                  if A[i] == x {
3                           swap(A[i], A[n-1]);
4                           break;
5                  }
6         return PARTITION(A, p, r);
}

在修改划分算法后，我们通过以下步骤来实现在n个元素的数组中找第i小元素的SELECT：
１、把数组A分成ceiling(n/5)个组，除最后一组外，每组都有5个元素，最后一组有n mod 5个元素；
２、对每组（的五个元素）用插入法进行排序，然后选出该组的中位数，即排好序的第三个元素；
３、对于步骤2中得到的所有的中位数，通过递归调用SELECT来找到它们的（下）中位数x，（也就是找到第2步得到的所有中位数中第floor(ceiling(n/5) / 2)小个元素）；
４、利用修改后的划分算法把元素分成小于x和大于x的两个子数组。如果设k为划分低区的元素个数加一，则x就是A中第k小的元素；
５、如果i = k，那我们就返回x，它便是我们要找的值。如果i < k，我们就在第4步里的划分低区继续递归调用SELECT来找到第i小的元素；如果i > k，我们就在划分高区递归调用SELECT找第i-k小的数。

下面是伪代码：

SELECT(A, p, r, i) {
// 步骤１、２
1      count = ceiling(n/5);
2         for i = 1 to count-1
3                  insertion sort A[(i-1)*5+1...i*5+1];
4         insertion sort A[(count-1)*5+1...n];
5         if count ==1
6                  return A[floor(n/2)];
// 步骤３
7         create array B;
9         for i = 1 to count-1
10  B[i] = A[(i-1)*5+3];
11  B[count] = A[(count-1)*5 + floor((n - (count-1)*5)/2)];
12  x = SELECT(B, 1, count, floor(count/2));
// 步骤４
13  q = PARTITION_X(A, p, r, i);
14    k = q-p+1;
// 步骤５
15  if i == k
16    return x;
17    elseif i < k
18    return SELECT(A, p, q-1, i);
19  else
20    return SELECT(A, q+1, r, i-k);
}

（上面的伪代码原书没有，博主根据书中描述得到的。）

下面分析下SELECT算法的性能。第2~4行中，共执行了ceiling(n/5)次排序操作，每次插入排序的代价都为O(1)，所以这几行总的代价为O(n)；第9~11行的创建数组B的代价为O(n)；第12行里，SELECT的输入规模为ceiling(n/5)，所以运行时间为T(ceiling(n/5))；第13行划分数组的代价为O(n)； 第18行，因为x是所有组里中位数的中位数，这样除最后组与包括x本身的组，至少有一半的组里有3个元素大于x，为3*(1/2 * floor(n/5) - 2) ≥ 3*n/10 -6，这样也意味着最多有n - (3*n/10-6) = 7*n/10 + 6个元素小于x，所以第18行低区元素的规模不会超过7*n/10+6；同理，第20行SELECT的输入规模也不会超过7*n/10+6。

综上所述，可以得到递归式：
T(n) ≤ T(floor(n/5)) + T(7*n/10+6) + O(n)

用代换法来解上面的递归式，设c, a为正常数得：
T(n) ≤ c*(floor(n/5)) + c*(7*n/10+6) + a*n
    ≤ c*(n/5) + c*(7*n/10+6) + a*n
    = 9*c*n/10 + 7*c + a*n
    = c*n + (-c*n/10 + 7*c + a*n)

为使T(n) ≤ c*n成立，-c*n/10+7*c+a*n ≤ 0必须成立。当n>70时，只需c ≥ 10*a*(n/(n-70)) 不等式成立。假设n > 140，有n/(n-70) ≤ 2，所以选c ≥ 20*a就可以。所以SELECT算法在最坏情况下运行时间为O(n)。