ソートと統計

ソートの話をするときは,はじめに軽くネタを振ってオーディエンスの偏り具合とかバラつき具合を調べておいて,最良の話題を選びましょう

Rの発言:

「ソートする時にさ、はじめにいっぺんデータを頭からなめて偏り具合とかバラつき具合を調べておいて、最良のソートアルゴリズムを選んで...」

これに周りのやつらが噛みつく:

「そんなもん、データ読み込むたんびにソートしとけば全部読み終わったときゃソート終わってんぢゃん」

R少なからず凹んだ様子。

まあデータを一回スキャンしたところで所詮 N Log(N) が N (Log(N)+1) になる程度の話なので (つまりマクロには無視できる),N が大きくなってくると普通に有効な戦略だと思います.別に全部調べなくてもサンプリングしても良いですし,並列化して事前スキャンという手もあるでしょう.
とはいえ,そういう話もオーディエンスは選んだ方が良いんでしょうな.N が小さい方面だと別の芸術と別の話題があるわけで.

quicksort と pivot 選択

quicksort の性能というか安定性が pivot の選び方の影響を受けるというのは,教科書でも割と強調されている話です.
例えば Haskell でいくら qsort が簡単に書けるといっても,毎回最初の要素を pivot に選んで本当に大丈夫なんでしょうか? と.

qsort []     = []
qsort (p:xs) = qsort [x | x <- xs, x < p] ++ [p] ++ qsort [x | x <- xs, x >= p]

んでまあ,quicksort を書きましょうという課題がさっさと終わって暇になった学生が次にやるのが,3 点ぐらいサンプリングしてきて真ん中の値をとるとか,もっとたくさんサンプリングしてきて真ん中の値を選ぶとかってのは実にありがち.これなんかは,パラメータチューニングとはいえ,「偏り具合とかバラつき具合を調べておいて,最良のソートアルゴリズムを選んで」の一種と言えるのではないでしょーか.
さて,対象の数列を均等に分割できる pivot を選べばパフォーマンスが良くなりそうだと直観的に思い付くところまでは良いとして,問題は「N 点サンプリング」した結果をどう使えばより均等に分割できるか,てとこですかね.
結局ここにいたって,対象の分布が分からないと何とも言えない問題にぶち当たるような気がします.一様分布なら N 点の平均を使えばうまくいきそうな気がしますが,一般論としてはどうでしょう?
突き詰めれば累積分布関数 (CDF) を考えて,CDF が 1/2 ぐらいのところをえいやと選ぶのが良さそうです.

分布関数と区間から計算で分割位置を求める

さて,そもそも最初から統計分布が分かっているのであれば,数列の min, max から適切な 分割基準が自動的に決まる気がしてきます.例えば相手が一様分布と分かっていれば,mix と max の平均を選べば良さそうです.LINQ の練習がてらにやってみましょう.

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

static class Program
{
    static IEnumerable<T> Repeat<T>(T value)
    {
        while (true) yield return value;
    }
    static void Main(string[] args)
    {
        // [0, 1] を返す疑似一様乱数ジェネレータ 
        var uniformGenerator = Repeat(new Random()).Select(r => r.NextDouble());

        // ジェネレータから一定数取り出す
        const int count = 100000;
        var seq = uniformGenerator.Take(count).ToArray();

        var min = seq.Min();
        var max = seq.Max();
        var threshold = (min + max)/2.0;
        var lessCount = seq.Count(x => x < threshold);
        Console.WriteLine("1st segment: {0}, 2nd segment: {1}", lessCount, count - lessCount);
    }
}

結果

1st segment: 49979, 2nd segment: 50021

当たり前のことではありますが,うまくいっているようです.並列計算でいきなり M 個のプロセッサに分散処理させたければ,区間を単純に M 分割すればよいでしょう.

分布関数と区間から計算で分割位置を求める (2)

一様分布では面白くないので,ロングテールで有名な Pareto 分布*1で実験してみましょう.
Wikipedia の 『Generating a random sample from Pareto distribution』を参考に一様分布から Pareto 分布を作ってみます.

static class Program
{
    static IEnumerable<T> Repeat<T>(T value)
    {
        while (true) yield return value;
    }
    static void Main(string[] args)
    {
        // [0, 1] を返す疑似一様乱数ジェネレータ 
        var uniformGenerator = Repeat(new Random()).Select(r => r.NextDouble());

        // K と MinX でパラメータ化された Pareto Distribution
        const double K = 2.0;
        const double MinX = 1.0;
        var paretoGenerator = uniformGenerator.Select(x => MinX * Math.Pow(x, -K));

        // ジェネレータから一定数取り出す
        const int count = 100000;
        var seq = paretoGenerator.Take(count).ToArray();

        var stat = new { Min = seq.Min(), Max = seq.Max() };

さて,ジェネレータから適当に何個か疑似乱数列を取り出して,その最大値が stat.Max, 最小値が stat.Min だったとします.今回はさすがにその平均で均等分割,とはいかなさそうです.
ここで Pareto 分布の累積分布関数に戻って計算しても良いですし,上の uniformGenerator の式の逆を考えても良いのですが,実は得られた数列の各項 x について Pow(x, -1/k) を計算すると,一様分布に戻ります.これを利用して,区間の均等分割が使えるというわけです.
簡単のために,Pow(stat.Min, -1/k) が 0 になり,Pow(stat.Max, -1/k) が 1 になるようにスケールを取り直してみましょう.

        // TODO: 最適化したい人は定数計算をループの外にくくりだして下さい
        var rescaledSeq = from x in seq
                          let newMin = Math.Pow(stat.Min, -1.0 / K)
                          let newMax = Math.Pow(stat.Max, -1.0 / K)
                          let newX = Math.Pow(x, -1.0 / K)
                          select (newX - newMin) / (newMax - newMin);

この rescaledSeq は [0, 1] に分布する一様分布になるわけですが,今ひとつぴんと来ない場合は uniformGenerator の定義と併せて考えてみてください.結局やっていることは単なる rand.NextDouble() になっているはずです.
再スケーリング結果は [0, 1] に分布しているので,適当に整数値 M をかけた後 int にキャストし,それをキーにグルーピングしてみましょう.これでだいたい要素数が揃った M 個の区間に分割できるはずです.

        const int NumSegment = 8;
        var groupedSeq = from x in seq
                         let newMin = Math.Pow(stat.Min, -1.0 / K)
                         let newMax = Math.Pow(stat.Max, -1.0 / K)
                         let newX = Math.Pow(x, -1.0 / K)
                         let rescaledX = (newX - newMin) / (newMax - newMin)
                         let segmentNumber = Math.Min((int)(NumSegment * rescaledX), NumSegment - 1)
                         group x by segmentNumber;

        foreach (var item in groupedSeq.OrderBy(pair => pair.Key))
        {
            Console.WriteLine("segment {0}: {1}", item.Key, item.Count());
        }

例えば 8 分割してみましょう.結果はこうなりました.

segment 0: 12499
segment 1: 12537
segment 2: 12448
segment 3: 12489
segment 4: 12648
segment 5: 12318
segment 6: 12424
segment 7: 12637

このように,形状母数が既知の Pareto 分布を分割統治でソートするときに,全スキャンまたはランダムサンプリングで区間の大きさを調べておくのは比較的有効と言えるんじゃないでしょうか.

まとめ

もちろん上の話は出来レース (ノイズの無い,形状母数が既知等) ではありますが,対象の統計分布が分かっていると色々と戦略の幅が増すのは事実です.特に今後は並列化をにらんで安定した処理分割が重要になってくると予想され,そこで対象の統計分布を利用しない手はありません.

そもそも人間活動や自然現象を相手にしていて,一様分布が出てくる方がむしろレア.例えば何かのデータを文字列長でソートするようなシナリオで,調べてみると 10 文字以下といった限られた範囲に 99% 以上の要素があったりするわけです.そんな状況でも,あなたは Array.Sort や std::sort を使いますか?
やはり今後は LINQ (to SQL) のように,what だけ書いておいて計算戦略は実行時に決定,という方向に期待したいところです.もちろん N が大きい世界の話.

*1:[http://reference.wolfram.com/mathematica/ref/ParetoDistribution.html:title=]