【Algorithm】选择问题与BFPRT算法

问题描述.

• 【选择问题】Selection Problem是求一个$n$个数的列表中第$k$小(实际上小还是大并不关键)的元素的问题，这个元素被称为第$k$个顺序统计量。
• 对于$k=1$或$k=n$的特殊情况，我们可以在线性时间内扫描该列表，然后从中选出最小(最大的元素)。
• 这个问题更加有趣的情况在于令$k=\lceil n/2\rceil$，该元素大于列表中半的元素，又小于其中另一半的元素，它被称为中值Median.
• 一个很直观的想法是，使用某种高效的排序算法对列表进行预处理，而后直接从中获得索引号为$k$的元素(假定从1开始计数)，使用诸如归并排序MergeSort或堆排序HeapSort这样高效又稳定的排序方法时，这一算法的时间复杂度为$O(n\cdot logn)+O(1)=O(n\cdot logn)$，不失为一种较好的方法。
• 但我们从信息量的角度来看，问题所要求的只是第$k$个最小的元素，实际上并不需要收集关于整个列表中$n$个元素的信息，所以$O(n\cdot logn)$很有可能不是这一问题解法所能达到的最好时间。

快速选择QuickSelect.

• 这一算法的思想类似于排序算法中的快速排序QuickSort，但也略有不同，也是选择一个枢轴量Pivot，使用某种划分算法将整个列表以其为基准划分，左边的都小于等于Pivot，右边的都大于等于Pivot，而后对于Pivot的位置进行判断：
  ①如果$Index[Pivot]=k$，那么Pivot就是第$k$顺序统计量；
  ②如果$Index[Pivot]>k$，那么就在Pivot的左边部分递归求解第$k$顺序统计量；
  ③如果$Index[Pivot<k]$，那么就在Pivot的右边部分递归求解第$k-Index[Pivot]$顺序统计量(因为从Pivot向左已经确定了前$Index[Pivot]$个元素)。
• 可以察觉到这一算法和快速排序的极高相似度，但某种程度上它又很类似于二分查找BinarySearch。划分算法的使用和快速排序很相似；而通过比较迭代结果与目标量来减治问题的规模这一点又和二分查找很类似。
• 我们暂且先不讨论具体的划分算法，而先关于快速选择算法的效率。对于一个n元素列表的划分，总是需要进行$n-1$次的键值比较，所以其最好情况下的效率就是$O(n)$；而如果好巧不巧，每一次迭代算法都产生一个<$0,n-1$>的划分组合，最坏情况下的效率就是${\sum }_1^{n-1}i\in O(n^2).$
• 一个相当复杂的算法分析表明，这种基于划分进行第$k$顺序统计量求解算法的平均效率是$O(n)$级别的。

划分算法.

• 【Hoare划分】Hoare是快速排序的发明者，英国计算机科学家，图灵奖获得者。
  
• 【Lomuto划分】这个人似乎不太出名，没有太多记载。
  

有关枢轴量的选择.

• 由于快速选择是一个基于划分的算法，所以有关快速排序中枢轴量选择的讨论也适合于快速选择。实际上后面要说的BFPRT算法就是基于稳定地选择一个好的枢轴量而提出的。
• 【随机枢轴量选择】即随机选择一个元素作为枢轴量，避免出现列表有序时<$0,n-1$>划分的尴尬情况。
• 【三平均划分】以数组最左边、最右边以及中间元素的中位数作为枢轴量，避免了点儿太背正好选到最小、最大元素的尴尬情况。

BFPRT算法.

• 上面给出的快速选择QuickSelect算法，能够在平均情况下达到线性时间复杂度，但我们还不够满意，我们希望找到一种最差时间复杂度更优的算法。
• 我们的确应该期待这样一个算法的存在，因为在第一部分给出的预排序再查找算法中，其最坏情况下的复杂度也不过是 $O(n\cdot logn)$，而选择问题中需要的第$k$顺序统计量需要的信息比排序算法还要少，所以我们有理由期待一个最坏情况下达到线性时间复杂度的算法。
• 幸运的是，这样的算法确实存在——BFPRT算法，这样一个略显怪异的名称其实是提出它的计算机科学家的名称首字母缩写，Blum、Floyd、Pratt、Rivest、Tarjan. 1973 年，五位大佬集体出动，合写了一篇题为《Time bounds for selection》的论文，在其中提到了线性第$k$顺序统计量查找算法，也称为中位数之中位数算法，而该算法的核心思想就在于优化枢轴量的选取，在每一次的迭代中都会丢弃Discard很大比例的元素，避免遇到那些导致 $O(n^2)$ 级别的坏划分出现。

算法详述.

• 为了不引起我个人翻译而导致的不准确，这里直接给出论文中该算法的描述(文中它被称为 $PICK$)，并且稍后会给出网络上几个关于BFPRT算法的描述版本，进行解释。

• 【说明】由于这段算法的叙述是有上下文的，所以其中有些符号显得有些没头没脑，以下是当中符号的说明(时间充裕的话还是推荐查看原论文)：
  
  
• 【四六级未过专属】$i$ $\theta$ $S$ 代表选取运算，即从集合 $S$ 中选取第 $i$ 小的元素；$x$ $\rho$ $S$ 代表查索引运算，即计算元素 $x$ 在集合 $S$ 中的次序；其中的变量 $b,c,d$ 都只是进行辅助作用，而网络上的几种不同说法的不一致就出现在这些辅助变量上。
• 步骤 $1(b)$ 中在选取枢轴量时，递归调用了 $PICK$ 算法，这样做的条件是递归调用该算法时不能只有一个分组，否则列表 $T$ 中仅有一个元素，何来的 $b$ $\theta$ $T$.

时间复杂度分析.

• 在第一步的分组操作中，需要对分好的每一个 $Column$ 中 $c$ 个元素进行排序，网络上关于这一步操作的时间复杂度说法都模模糊糊，要么就是一笔带过，要么就是压根不提，最后也给出一个BFPRT算法能够在最坏情况下达到 $O(n)$ 的结论。以我的认知，比较排序算法的渐进复杂度下界是 $O(n\cdot logn)$，某些非比较排序在适合的情况下会达到线性级别，但整个BFPRT算法的最差复杂度被证明为 $O(n)$，这其中必然存在值得深究的东西。
• 文中的论述如下所示：
  
  这里出现的代价 $h(c)$ 并不是渐进情况下假定规模趋近于无穷的，当中提到两篇论文，一篇关于 $Ford$ - $Johnson$ 算法(我并没有找到这两篇论文，可能是年代久远了)，另一篇中给出了比较 $c$ 个数的代价，也就是上图中的那个公式(至于为什么，反正我现在是没搞懂)。所以操作 $1(a)$ 的代价是 $n\cdot h(c)/c$ ，使得 $c$ 必须是一个以常量为上界的辅助量，否则即使 $c$ 是一个线性量，$h(c)$ 所代表的求和式近似为 $log(n!)$ ，也就近似为 $n\cdot logn$.
• 记$P(n)$为算法 $PICK$ 的最差情况下时间复杂度，在操作 $1(b)$ 中我们有这样一个运算 $m=b$ $\theta$ $T$ ，其中 $T$ 是 $n/c$ 个 $Column$ 中第 $d$ 小元素的集合，所以 $P(n)$ 中至少包含了 $P(n/c)$.
• 选择好枢轴量 $m$ 之后，整个列表的情况可以通过下图来理解：
  
• 以$m$为基准，这 $n$ 个数被分成了四个 $Block$ ，很容易理解，$G$ 的块中元素全部大于 $m$，而 $L$ 的块中元素全部小于 $m$，所以在操作 $Ste{p}_2$ 中需要和 $m$ 进行比较的元素只有块 $A$ 和块 $B$，就能够确定出 $m$ $\rho$ $S$. 在最坏的情况下，$m$ 所面临的待比较的元素个数不过是整个列表的元素个数 $n$，所以本步骤的代价为 $n$ .
• 操作 $Ste{p}_3$ 中所做的是根据 $Ste{p}_2$ 中计算出的 $m$ $\rho$ $S$ 来决定下一次迭代时的参数。如果 $m$ 的位置就是那个指定位置，那么算法结束，输出所有比 $m$ 小的数以及 $m$ 本身；如果 $m$ 的位置大于指定位置，说明 $m$ 过大，递归地对于除去大于等于 $m$ 所有数后的集合调用算法；否则递归地对于除去小于等于 $m$ 所有数后的集合调用算法，注意此时的参数 $i$ 需要改变。这一步骤中，最差的情况就是仅仅除去了块 $G$ 和 $L$ 中较小的那一个，而后对于剩余集合递归调用算法，这已经是最差的划分情况了，所以 $Ste{p}_3$ 的代价为 $P(n-min(|G|,|L|))$ .
• 综上所述，整个 $PICK$ 的代价，也就是时间复杂度有下式：
• $P(n)\le n\cdot h(c)/c+P(n/c)+n+P(n-min(|G|,|L|))$.
• 注意这个表达式中依旧是存在辅助量的，实际使用算法时我们需要给出这些辅助量。论文中使用的一组辅助量值是($c:21,d:11,b:n/2c$)，所以上式可以写为：
• $P(n)\le n\cdot h(21)/21+P(n/21)+n+P(31n/42)$.
• 由于 $h(21)=66$，最后通过数学方法解得 $P(n)\le 58n/3\in O(n)$，确实是一个最坏情况为 $O(n)$ 的算法。另外，我们需要确保 $c\ge 5$，否则会导致整个形为 $c\times (n/c)$ 的矩阵过于扁平，每一次迭代前无法除去足够多的元素而失去线性的时间复杂度。
• 【总结】宏观上来说，BFPRT算法和快速选择算法很类似，都是在划分整个列表后，根据枢轴量的位置进行筛选、决策，而后开始下一步的迭代。

《Time Bounds for Selection》： PICK is quite similar to the algorithm FIND (Hoare [5]), except that the element m about which to partition S is chosen more carefully.

• 关于当中提到的Ford-Johnson算法，可以参考这篇论文中的描述，个人觉得这个算法就是锦标赛排序算法。

来自网络.

• 摘自知乎。
  
• 这一BFPRT算法的描述版本是出现较多的，其中指定了辅助量 $c=5$，关于步骤2和步骤3中的中位数选取，实际上就是指定辅助量 $d=(c+1)/2$ 以及 $b=n/2c$。但步骤3中的最后一个条件【直到剩下一个数字】这一描述不是很直观，他想要表达意思应该是递归调用到某一次执行算法时第2步只得到了一个中位数就停止。如果表述为【递归调用，直到 $n/5\le 1$ 会更清楚】。
• 后面也有人指出了这一问题，发言摘录如下，暴躁老哥，在线怼人：
  

涉及算法のC++实现.

快速排序QuickSort.

#include <iostream>
#include <vector>
#include <algorithm>
#define Limit 5

using namespace std;

typedef int ElementType;

typedef vector<ElementType> Line;

void QuickSort(Line& line,int left,int right);

void QuickSortSub(Line& line, int left, int right);

int Partition(Line& line, int left, int right);

void Display(Line& line);

int main()
{
	Line line;
	int number = 0;
	while (cin >> number)
	{
		if (cin.get() == '\n')
		{
			break;
		}
		line.push_back(number);
	}
	line.push_back(number);

	QuickSort(line, 0, line.size());

	Display(line);
	
	return 0;
}

void QuickSort(Line& line,int left,int right)
{
	if (right - left <= 2) //Too less elements.
	{
		//Use some simple sort algorithm.
		//I am lazy,so i call C++.sort().
		sort(line.begin(), line.end());
	}
	QuickSortSub(line, left, right);
}

void QuickSortSub(Line& line, int left, int right)
{
	if (right - left <= 2)
	{
		sort(line.begin(), line.end());
		return;
	}
	int partition = Partition(line, left, right);
	QuickSortSub(line, left, partition);
	QuickSortSub(line, partition + 1, right);
}

int Partition(Line& line, int left, int right)
{
	int pivot = line[left];
	int temp = 0;

	int i = left + 1;
	int j = right - 1;

	while (i < j)
	{
		while (line[i] < pivot)
		{
			++i;
		}

		while (line[j] > pivot)
		{
			--j;
		}
		temp = line[i];
		line[i] = line[j];
		line[j] = temp;
	}
	temp = line[j];
	line[j] = line[i];
	line[i] = temp;

	line[left] = line[j];
	line[j] = pivot;

	return j;
}

void Display(Line & line)
{
	for (int i = 0; i < line.size(); ++i)
	{
		cout << line[i] << "    ";
	}
}

归并排序MergeSort.

堆排序HeapSort.

快速选择QuickSelect.

PICK算法BFPRT.

• $PICK$这个名字是五位计算机科学家在论文中的称呼。

#include <iostream>
#include <vector>
#include <algorithm>
#define Limit 5

using namespace std;

class Tuple
{
public:
	int key;
	int value;
	Tuple(int a, int b) :key(a), value(b)
	{

	}
};

typedef int ElementType;

typedef vector<ElementType> Line;

typedef vector<Tuple> BinaTuple;

void BFPRT(Line line,int total,int index,int left,int right);

int EfficientPartition(Line& line,int total,int left,int right);

//Compare function.
bool myfunction(Tuple i, Tuple j) 
{ 
	return (i.value < j.value);
}

int main()
{
	Line line;
	int number = 0;
	while (cin >> number)
	{
		if (cin.get() == '\n')		//If EOF,end.
		{
			break;
		}
		line.push_back(number);
	}
	line.push_back(number);
	int index = 0;
	cin >> index;
	BFPRT(line,line.size(),index,0,line.size());
	return 0;
}

void BFPRT(Line line,int total,int index,int left,int right)
{
	int partition = EfficientPartition(line, total, left, right);
	//Get the position of pivot.

	//Decrease-and-conquer by compare partition and index.

	if (partition == index - 1)		//Equal:print the value.
	{
		for (int i = 0; i < index; ++i)
		{
			cout << line[i] << "    ";
		}
	}
	else if (partition > index - 1) //Partition too large:discard elements too large.
	{
		line.erase(line.begin() + partition, line.begin() + total);
		BFPRT(line, partition, index, 0, partition);
	}
	else		//Partition too small:discard elements too small,and index should be reseted.
	{
		for (int i = 0; i < partition + 1; ++i)
		{
			cout << line[i] << "    ";
		}
		line.erase(line.begin(), line.begin() + partition + 1);
		BFPRT(line, total - partition - 1, index - partition - 1, 0, total - partition - 1);
	}
}

int EfficientPartition(Line& line, int total, int left, int right)
{
	//Scale of elements too small,dispose directly.
	if (total == 2 || total ==1)
	{
		sort(line.begin(), line.end());
		return 0;
	}

	//Divide into groups.
	int group = total / Limit;
	BinaTuple temp;

	for (int i = 0; i < group; ++i)
	{
		sort(line.begin() + i * Limit, line.begin() + (i + 1) * Limit);
		temp.push_back(Tuple(i*Limit + Limit / 2, line[i*Limit + Limit / 2]));
		//Choose median in each group.
	}

	//The last group.
	if (total%Limit != 0)
	{
		sort(line.begin() + group * Limit, line.end());
		int mid = (group*Limit + total - 1) / 2;
		temp.push_back(Tuple(mid, line[mid]));
	}

	sort(temp.begin(), temp.end(), myfunction);
	
	//Get the median of median:pivot.
	int mid = temp.size() / 2;
	int pivot = temp[mid].value;
	int index = temp[mid].key;

	//Put pivot at the front,SWAP(line[0],pivot).
	int hoe = line[left];
	line[left] = pivot;
	line[index] = hoe;

	//Compute the eventual position of pivot.
	//HoarePartition.
	int i = left + 1;
	int j = right-1;
	while (i < j)
	{
		while (line[i] < pivot)
		{
			++i;
			if (i == right)
			{
				--i;
				break;
			}
		}
		while (line[j] > pivot)
		{
			--j;
			if (j == left)
			{
				++j;
				break;
			}
		}
		//SWAP(line[i],line[j]).
		hoe = line[i];
		line[i] = line[j];
		line[j] = hoe;
	}
	//Rollback the last SWAP().
	hoe = line[j];
	line[j] = line[i];
	line[i] = hoe;

	//SWAP(line[j],pivot)
	line[0]=line[j];
	line[j]=pivot;
	
	//j is the eventual position of pivot.
	return j;
}