· 開発

中央値の中央値（BFPRT）：最悪の場合でも線形時間でk番目の値を求める

1973年のBlum・Floyd・Pratt・Rivest・Tarjanの論文にあるPICKアルゴリズムが最悪の場合でも線形時間でk番目の値を選べる理由と、5個ずつのグループ分けの由来を原論文に沿って整理します。

algorithmsselectionmedian-of-medianscomplexity

n個の値のうちk番目に小さい値を求める問題を選択（selection）問題と呼びます。代表的な例は中央値です。全体をソートしてk番目を取り出せば答えは出ますが、値を一つ得るために全要素の順序を決めていることになります。

Blum、Floyd、Pratt、Rivest、Tarjanは1973年の論文「Time Bounds for Selection」で、kが何であっても比較回数がnの線形関数を超えないアルゴリズムPICKを示しました。ピボットの選び方から一般に「中央値の中央値（median of medians）」と呼ばれ、著者の頭文字からBFPRTとも呼ばれます。この記事では、この論文（Journal of Computer and System Sciences 第7巻 448–461ページ、共著者Rivestが公開している写し）に基づいて、PICKの構造、線形になる理由、そして「なぜ5個ずつ分けるのか」という疑問を扱います。その後の改良研究や特定の言語ライブラリの実装は扱いません。

ピボット一つで性能が決まる

k番目の値を求めるいちばん素直な方法は、クイックソートの分割を応用することです。ピボットmを一つ選んで残りの値とすべて比較し、k番目の値があり得ない側を捨てて、残った側で繰り返します。論文は、PICKはHoareのFINDアルゴリズムとよく似ており、分割の基準となるmをより慎重に選ぶ点だけが違うと説明しています。

この方法の弱点はピボットです。以下はこの記事による計算です。ピボットが毎回、残りの値の最小値か最大値だと、一回の分割で値が一つずつしか減りません。すると比較回数はおよそ n + (n − 1) + … となり、n²/2 程度まで増えます。クイックソートでピボットの選び方が最悪の場合を生む過程は、クイックソートの最悪の場合の記事で扱いました。PICKはピボット選びにより多くの比較を使い、分割のたびに一定の割合以上を捨てることを保証します。

PICKの手順

論文はPICKを三つのパラメータ c、d、b で一般化して説明しています。集合Sからi番目に小さい値を求めるとき、1回の反復は次のとおりです。

  1. Sを長さcの列 n/c 本に分け、各列をソートします。各列のd番目に小さい値を集めた集合をTとし、Tのb番目の値をPICKの再帰呼び出しで求めてmとします。
  2. mとの大小関係がまだわからない要素をすべてmと比較し、S内でのmの順位を求めます。
  3. 順位がiならmが答えなので停止します。順位がiより大きければm以上の要素を、小さければm以下の要素を捨て（このときiも捨てた個数だけ減らします）、手順1に戻ります。

dを列の真ん中、bをTの真ん中に取ると、mは「各列の中央値の中央値」になります。名前はここから来ています。論文は、こうして捨てる部分集合が毎回、残りの要素の4分の1以上を含むと述べています。

線形になる理由

論文は最大比較回数をP(n)とし、c = 21、d = 11、b = n/42 を選んで次の漸化式を立てます（式6）。

P(n) ≤ 66n/21 + P(n/21) + n + P(31n/42)

各項の意味は次のとおりです。66はFord–Johnsonの方法で21個をソートする比較回数h(21)なので、第1項は n/21 本の列をすべてソートするコストです。P(n/21)は n/21 個の列中央値から中央値を再帰で求めるコスト、nはmと残りの要素を比較するコストの上限、最後のP(31n/42)は捨てた後に残った要素に対する次の反復です。論文はここから帰納法で P(n) ≤ 58n/3 を導きます。

31n/42がどこから来るかは、論文の図1から読み取れます。以下はその図に基づくこの記事の計算で、列がちょうど割り切れない場合や±1の差は無視しています（論文も、nがcで割り切れない場合の細部は省略すると明記しています）。列中央値の半分、つまり n/42 本の列は中央値がm以下です。列はソート済みなので、これらの列では下側の11個がすべてm以下であり、m以下だと確実にわかる要素が11n/42個以上あります。m以上の側も同じ理由で11n/42個以上です。どちらの側を捨てても11n/42個以上が消えるので、残る要素は多くても31n/42個です。

線形時間になる鍵は、二つの再帰呼び出しのサイズの和です。n/21 + 31n/42 = 33n/42 で、nより小さくなります。次の段階で扱う要素数の合計が毎回一定の比率で減るので、全段階のコストを足すと収束する等比級数になり、全体がnに比例します。この部分はこの記事による説明で、論文は同じ内容を「各段階でSのかなりの割合を捨て、そのコストが捨てた要素数に比例するから」とまとめています。

5個ずつという数はどこから来るのか

中央値の中央値を「5個ずつに分ける」と習ったなら、原論文の数字は少し意外かもしれません。基本の証明は21個の列を使い、定数を小さくした改良版PICK1は15個の列を使います。5という数は、証明の最後の一文にしか出てきません。「PICKが線形時間で動くには c ≥ 5 でなければならないことに注意せよ。」論文はこの一文に別途の証明を付けていません。

理由は、前の節の計算をcについて一般化すると見えてきます。これもこの記事による計算です。真ん中の要素が一つに決まるようにcを奇数 2k + 1 とすると、m以下だと確実にわかる要素は列の半分からそれぞれ k + 1 個ずつなので約 (k + 1)n/(2c) 個で、残る要素は多くても n − (k + 1)n/(2c) 個です。二つの再帰呼び出しのサイズの和をnに対する比率で書くと次のようになります。

列のサイズ c 中央値の再帰 次の反復 和
3 1/3 2/3 1
5 1/5 7/10 9/10
7 1/7 5/7 6/7
21 1/21 31/42 11/14

c = 3 では和がちょうど1なので、段階を下っても扱う要素数の合計が減らず、前の節の等比級数の議論では線形の上限が得られません。5は和が1未満になる最小の奇数で、5個ずつ分ける説明はこの最小値を使っていることになります。ただし、この表はPICKの構造をそのまま保った単純な上限の計算であり、構造を変えた別の変形について判定するものではありません。

逆に、cをいくらでも大きくできるわけでもありません。手順1の列のソートにかかるコストは n·h(c)/c で、論文はPICKが線形であるためにはcが定数で抑えられていなければならないと書いています。論文が基本の証明でc = 21を選んだのも、漸化式(5)の上限を最小にするためでした。同じ漸化式にc = 5を入れると（5個のソートに比較7回、次の反復は7n/10）上限は24nとなり、c = 21のときの58n/3より大きくなります。この比較はこの記事による計算です。5は成り立つ最小の値であって、定数を小さくするのに最適な値ではありません。

論文が示した定数と限界

論文の結論は、選択は最悪の場合でも線形時間で行えるということです。改良版PICK1では、どのiについても比較回数が約5.43nを超えないことを証明しており（正確には5.4305nで、最後の5が循環）、iが両端に近いときはPICK2でさらによい上限を得ています。下限については、中央値を求めるには少なくとも 3n/2 − 2 回の比較が必要であることを示しました。

実用面については、論文自身が但し書きを付けています。PICK1は付随的な処理が多すぎて（原文は “full of red tape”）実際に効率よく実装することはできないと書き、上限と下限のどちらの定数も大きく改善できるだろうと見ていました。論文の分析対象は比較回数で、入力は互いに異なるn個の数と仮定しています。

Pythonに移した例

以下のコードは、5個ずつ分けるPICKをPythonに移したものです。論文は互いに異なる値を仮定していますが、このコードはピボットと等しい値を別に数えることで重複する値も扱えるようにしています。kは0から数えます。

def select(a, k):
    """a の中で k 番目（0 から数える）に小さい値を返す。"""
    if len(a) <= 5:
        return sorted(a)[k]
    # 5 個ずつのグループに分け、各グループの中央値を集める（最後のグループは短いことがある）
    groups = [sorted(a[i:i + 5]) for i in range(0, len(a), 5)]
    medians = [g[len(g) // 2] for g in groups]
    # 中央値の中央値を再帰で求めてピボットにする
    pivot = select(medians, len(medians) // 2)
    lo = [x for x in a if x < pivot]
    hi = [x for x in a if x > pivot]
    n_eq = len(a) - len(lo) - len(hi)
    if k < len(lo):
        return select(lo, k)
    if k < len(lo) + n_eq:
        return pivot
    return select(hi, k - len(lo) - n_eq)


print(select([9, 1, 8, 2, 7, 3, 6, 4, 5, 0, 5], 5))  # 5

このコードは構造を見せるためのものです。段階ごとに新しいリストを作り、論文のように結果がすでにわかっている比較を省くこともしていません。比較回数の定数は論文の分析とは異なりますが、二つの再帰呼び出しのサイズの和がnより小さくなる構造は同じです。

確認しておきたいこと

  • 保証の種類：PICKの線形の上限は、確率や平均ではなく最悪の場合についての決定的な保証です。
  • グループのサイズ：原論文は5に固定していません。5はこの構造で線形の議論が成り立つ最小の奇数で、論文の基本の証明は21、改良版は15を使っています。
  • 定数：論文が証明した上限は約5.43n回の比較で、著者自身もこの定数には大きな改善の余地があると見ていました。
  • 重複する値：論文は互いに異なる値を仮定しています。実装ではピボットと等しい値を別に扱う必要があります。

ブログへ