C++で数値ストリームの中央値を効率的に求める方法

この問題では、整数が継続的に読み込まれていくデータストリームが与えられます。目的は、要素を読み込みながら、その時点までの要素群に対して中央値を計算していくプログラムを作成することです。

中央値とは

中央値(Median)とは、ソートされた数列(昇順・降順のいずれでも可)において中央に位置する要素のことです。

• 要素数が奇数の場合：中央値は中央の要素そのもの
• 要素数が偶数の場合：中央値は中央に位置する2つの要素の平均値

具体例で理解する

入力として「3, 65, 12, 20, 1」が与えられた場合、各入力の時点での中央値は次のように変化していきます。

Input - 3 : sequence -(3) : median - 3
Input - 65 : sequence -(3, 65) : median - 34
Input - 12 : sequence -(3, 12, 65) : median - 12
Input - 20 : sequence -(3, 12, 20, 65) : median - 16
Input - 1 : sequence -(1, 3, 12, 20, 65) : median - 12

ここでは説明を分かりやすくするため、ソート済みの数列を表示しています。

解決アプローチ

この問題には主に3つの解法があります。

1. 毎回ソートを実行してから中央値を求める方法
2. 自己平衡二分探索木(セルフバランシングBST)を利用する方法
3. ヒープを利用する方法

この中で最も有望なのがヒープを使う手法です。最大ヒープと最小ヒープを組み合わせることで、要素を挿入するたびに効率よく中央値を取得できます。

基本的な考え方は以下の通りです。

• これまで読み込んだ数のうち、小さい方の半分を最大ヒープ(左側)に格納する
• 大きい方の半分を最小ヒープ(右側)に格納する
• 中央値は常に、両ヒープの先頭要素(またはその平均)から求められる

この方法なら、挿入操作はO(log n)、中央値の取得はO(1)で行えるため、ストリーム処理に非常に適しています。

C++での実装例

上記の解法をC++で実装したプログラムがこちらです。

#include <iostream>
using namespace std;
#define MAX_HEAP_SIZE (128)
#define ARRAY_SIZE(a) sizeof(a)/sizeof(a[0])
void swap(int &a, int &b){
    int temp = a;
    a = b;
    b = temp;
}
bool Greater(int a, int b){
    return a > b;
}
bool Smaller(int a, int b){
    return a < b;
}
int Signum(int a, int b){
    if( a == b )
        return 0;
    return a < b ? -1 : 1;
}
class Heap{
    public:
        Heap(int *b, bool (*c)(int, int)) : A(b), comp(c)
        { heapSize = -1; }
        virtual bool Insert(int e) = 0;
        virtual int GetTop() = 0;
        virtual int ExtractTop() = 0;
        virtual int GetCount() = 0;
        protected:
        int left(int i) {
            return 2 * i + 1;
        }
        int right(int i) {
            return 2 * (i + 1);
        }
        int parent(int i) {
            if( i <= 0 ){
                return -1;
            }
            return (i - 1)/2;
        }
        int *A;
        bool (*comp)(int, int);
        int heapSize;
        int top(void){
            int max = -1;
            if( heapSize >= 0 )
                max = A[0];
            return max;
        }
        int count() {
            return heapSize + 1;
        }
        void heapify(int i){
            int p = parent(i);
            if( p >= 0 && comp(A[i], A[p]) ) {
                swap(A[i], A[p]);
                heapify(p);
            }
        }
        int deleteTop(){
            int del = -1;
            if( heapSize > -1){
                del = A[0];
                swap(A[0], A[heapSize]);
                heapSize--;
                heapify(parent(heapSize+1));
            }
            return del;
        }
        bool insertHelper(int key){
            bool ret = false;
            if( heapSize < MAX_HEAP_SIZE ){
                ret = true;
                heapSize++;
                A[heapSize] = key;
                heapify(heapSize);
            }
            return ret;
        }
};
class MaxHeap : public Heap {
    private:
    public:
        MaxHeap() : Heap(new int[MAX_HEAP_SIZE], &Greater) { }
        int GetTop() {
            return top();
        }
        int ExtractTop() {
            return deleteTop();
        }
        int GetCount() {
            return count();
        }
        bool Insert(int key) {
            return insertHelper(key);
        }
};
class MinHeap : public Heap{
    private:
    public:
        MinHeap() : Heap(new int[MAX_HEAP_SIZE], &Smaller) { }
        int GetTop() {
            return top();
        }
        int ExtractTop() {
            return deleteTop();
        }
        int GetCount() {
            return count();
        }
        bool Insert(int key) {
            return insertHelper(key);
        }
};
int findMedian(int e, int &median, Heap &left, Heap &right){
    switch(Signum(left.GetCount(), right.GetCount())){
        case 0: if( e < median ) {
            left.Insert(e);
            median = left.GetTop();
        }
        else{
            right.Insert(e);
            median = right.GetTop();
        }
        break;
        case 1: if( e < median ){
            right.Insert(left.ExtractTop());
            left.Insert(e);
        }
        else
            right.Insert(e);
            median = ((left.GetTop()+right.GetTop())/2);
        break;
        case -1: if( e < median )
            left.Insert(e);
        else {
            left.Insert(right.ExtractTop());
            right.Insert(e);
        }
        median = ((left.GetTop()+right.GetTop())/2);
        break;
    }
    return median;
}
void printMedianStream(int A[], int size){
    int median = 0;
    Heap *left = new MaxHeap();
    Heap *right = new MinHeap();
    for(int i = 0; i < size; i++) {
        median = findMedian(A[i], median, *left, *right);
        cout<<"Median of elements : ";
        for(int j = 0; j<=i;j++) cout<<A[j]<<" ";
        cout<<"is "<<median<<endl;
    }
}
int main(){
    int A[] = {12, 54, 9, 6, 1};
    int size = ARRAY_SIZE(A);
    printMedianStream(A, size);
    return 0;
}

出力結果

Median of elements : 12 is 12
Median of elements : 12 54 is 33
Median of elements : 12 54 9 is 12
Median of elements : 12 54 9 6 is 10
Median of elements : 12 54 9 6 1 is 9

まとめ

本記事では、連続的に流れてくる整数ストリームから中央値をリアルタイムに求める方法を紹介しました。毎回ソートし直す方法(O(n log n))と比べて、最大ヒープと最小ヒープを併用する手法は挿入ごとにO(log n)で済むため、大量のデータを扱うストリーム処理において大幅な性能向上が期待できます。