P99与蓄水池算法(reservoir sampling)
一个监控问题
假设你是一个运维工作人员，维护着一个访问量巨大的服务，然后有一天，老板跑来问你这个服务的 p99 响应时间是多少？（p99 响应时间：系统 99%的请求都快于这个时间，而 1%的请求则慢于这个时间。即响应时间的 99% 分位点。）
一个简单的算法是直接对所有请求的访问时间进行排序，然后取出 99%位置的请求的访问时间。但是由于系统访问量巨大，每天都有上亿的请求，而一般的排序算法的时间复杂度为 O(n * log(n))，在这样的数据量级上，算法的空间复杂度也很高。怎么优化这个算法呢？
抽样
上亿的数据实在太多了，但是假设我们的系统是相对稳定的，则所有数据应该是同分布的，那我们应该可以通过随机采样的方式，选出一部分样本来，用这些样本来描述总体。比如说我知道今天会有一亿个请求，然后我就可以按千分之一来进行数据抽样，即每个数据都按照千分之一的概率抽取到样本集中。
这里就遇到另一个关键性的问题，我不知道今天会有多少请求，即我不知道总体有多大，更讨厌的是，总体在不断变大，那样本集也会不断变大，排序会越来越困难。
这时候，你应该打开 John Bentley 所著的《编程珠玑》，翻到第十二章，那里有这个问题的答案。
答案就是蓄水池算法，具体可参考《编程珠玑》或 wiki。https://en.wikipedia.org/wiki/Reservoir_sampling , 算法如下。
for (long i = 0; i < LOOP_NUM; i++) {
if (i < SAMPLE_NUM) {
samples.push_back(n(e));
} else {
long r = randint(0, i);
float ne = n(e);
if (r < SAMPLE_NUM) {
samples[r] = ne;
}
}
}
注意，这个 samples 数组是动态更新的，每个新样本进来都会被随机选择是否替换进来。
排序
回到上面的 p99 问题，我们可以设定程序每隔一段时间对样本数组进行排序，并选择 99%位置的数字。时间复杂度就是 O(m * log(m))，空间复杂度为 O(m)，其中 m 为样本集大小。
但是，我们还有更好的方法。这显然是一个 Top k 问题，看到 Top k 就想到堆嘛。
我们可以维护一个容量为 k 的最小堆（这里 k = 0.1 * m），然后把元素一个一个加入堆，当堆满了之后，只有比堆顶元素大的元素才能加入堆，并顶掉原来的堆顶元素。最后在堆顶的元素就是 p99 的值啦。这样的方法在时间和空间上都能省不少。（另外，如果你想知道最快的 1%的话，就要用最大堆了。）
//最小堆
int q1_cap = (int) (SAMPLE_NUM * 0.01);
for (int i = 0; i < SAMPLE_NUM; i++) {
if (i < q1_cap) {
q1.push(v[i]);
} else {
if (v[i] > q1.top()) {
q1.push(v[i]);
q1.pop();
}
}
}
Put it together
我们做个实验，实验数据不是系统延时，而是一个正态分布样本生成器。我们用样本生成器生成一亿条数据，并按千分之一进行采样，然后利用堆计算其 p99（其实就是 1%分位点。）
#define LOOP_NUM 100000000
#define SAMPLE_NUM 100000
#define randint(min, max) (rand() % (max - min + 1) + min)
int main() {
priority_queue<float, std::vector<float>, std::greater<float>> q1; // 最小堆
std::normal_distribution<float> n(0, 1); //均值为 0, 方差为 1
std::default_random_engine e; //random engine
std::vector<float> v;
clock_t t1 = clock();
// reservoir sampling
for (long i = 0; i < LOOP_NUM; i++) {
if (i < SAMPLE_NUM) {
v.push_back(n(e));
} else {
long r = randint(0, i);
float ne = n(e);
if (r < SAMPLE_NUM) {
v[r] = ne;
}
}
}
clock_t t2 = clock();
//最小堆
int q1_cap = (int) (SAMPLE_NUM * 0.01);
for (int i = 0; i < SAMPLE_NUM; i++) {
if (i < q1_cap) {
q1.push(v[i]);
} else {
if (v[i] > q1.top()) {
q1.push(v[i]);
q1.pop();
}
}
}
clock_t t3 = clock();
std::cout << q1.top() << std::endl;
std::cout << "sampling time " << (t2 - t1) / 1000000. << std::endl;
std::cout << "heap time " << (t3 - t2) / 1000000. << std::endl;
return 0;
}
最终结果为 2.32081。堆化非常快，只用了 2.5ms。我们知道标准正态分布的上 0.01 分位点是 2.327，已经非常接近了。
同时请注意，我们的采样算法不需要知道总体样本数，也就是说，我们继续向其中加入一亿数据，应该还是非常接近 2.327 的。我们试试在上面的 return 前加入如下代码：
....
for (long i = 0; i < LOOP_NUM; i++) {
float ne = n(e);
if (i < SAMPLE_NUM) {
v.push_back(ne);
} else {
long r = randint(0, i + LOOP_NUM);
if (r < SAMPLE_NUM) {
v[r] = ne;
}
}
}
for (int i = 0; i < SAMPLE_NUM; i++) {
if (i < q1_cap) {
q1.push(v[i]);
} else {
if (v[i] > q1.top()) {
q1.push(v[i]);
q1.pop();
}
}
}
std::cout << q1.top() << std::endl;
return 0;
}
估算出来的 p99 为 2.32563，离 2.327 又近了一步。
PS: 代码可以在我的 github 找到，下面这个仓库的 sample 目录下 github.com/chenminhua/math