変なソートを作りたい

この記事は呉高専アドベントカレンダー2022に11月30日まで空きがあったら適当に埋めるための記事です。→ 4日目の記事になりました ~~/ 空きがないので放流します~~。

ごあいさつ
はじめに
既存の有名ソートアルゴリズムたち
変なソートを作ろう
- 平方分割バブルソート
  - 追記
- 最長増減部分列分割ソート
  - 最長増加部分列
実験
- ソースコード
おわりに

ごあいさつ

こんにちは、点Nです。電気情報工学科卒です。高専在学中は高専プロコンに出たり、呉高専史上初めてICPC（大学対抗競プロ合戦）に参加して予選敗退したりといった感じです。

こちらは11日の記事の息抜きに書いていたものです。さっそく始めます。

はじめに

ソート。配列を昇順（あるいは降順）に並び替える処理です。単純ながらあらゆる場面で頻繁に使用するアルゴリズムであり、その実現方法には様々なものがあります。呉高専OBである私も、情報処理の授業で各種ソートを実装した覚えがあります。（余談：当時はC言語でしたが、シラバスを見たところ今の子たちはPythonでこれらを学ぶっぽいです。またシラバスにはグラフ理論や動的計画法（ナップサック問題）、乱択アルゴリズムなどに関する記述もあり、時代の変化を感じます）しかし、教科書に書いてあることを写してソート関数完成やったーというのは（必要かもしれませんが）あまり楽しいものではありません。自分でソートのアルゴリズムを考えて実装したほうが楽しいはずです。なのでそれを今からやりたいと思います。まずアルゴリズムの教科書によく書いてある有名なソートをいくつか振り返ります（知っている人は飛ばしてください）。次に、変なソートを作ってみて、その性能について考えてみたいと思います。途中でちょっとだけ競技プログラミング（競プロ）の話もします。

お約束

ソートしたい配列をA、その長さをNとします。0-index（AはA[0]からA[N-1]まで）です。ソートは常に昇順です。
計算量はたいてい平均時間計算量の意味です。計算量の表現にオーダ記法を用います。使い方が雑ですが目をつぶってください。知らない人は、中に書かれた関数をプロットしたときに急に増えていかないほうが嬉しい、と思ってください。
Aの連続した区間は半開区間でA[a:b)のように表します。aを含みbを含まない、A[a]からA[b-1]まで、の意味です。

既存の有名ソートアルゴリズムたち

バブルソート

Aを走査しながら、A[i]>A[i+1]となっている部分をswapします。走査回数が $N$ 回、走査する部分の長さが平均 $N/2$ なので、計算量は $O(N ^ {2})$ です。単純かつ実装難易度が低いですが、計算量は（実用に耐えうる）ソートアルゴリズムの中では悪く、優れたアルゴリズムではありません。

選択ソート

A[0:x)がソートされているとします。未ソートの範囲A[x:N)を走査します。最小の値がA[min]にあるとして、これをA[x]とswapするとA[0:x+1)がソート済みになります。これを繰り返します。計算量はバブルソートと同じ理屈で $O(N ^ {2})$ です。

クイックソート

正直この記事とあまり関係ないので紹介しなくてもいいけどこれ紹介しないのは「アンパンマンワールドのゆかいななかまたちを紹介します！」って言ってアンパンマンを紹介しないのと同じようなものじゃんという感じなので仕方なく紹介します。

A[0:N)をソートすることを考えます。適当な値をピボットとして設定し、それ未満の値をA[0:x)、それ以上の値をA[x:N)の範囲に持ってきます。A[0:x)とA[x:N)がソートされれば全体をソートしたことになるので、これで小さな2つの問題に分割されました。この2つの問題も、同じようにそれぞれ小さな2つの問題に分割して解きます。問題の最小単位は長さ2以下の配列に対するソート（単純な大小比較とswap）です。

検索すればいくらでも出てくるので理屈は他の文献に委ねますが、計算量は平均 $O(N\log N)$ 、かつ同程度に速い他の手法と比較しても優れていると言われる優秀なアルゴリズムです。しかし、最悪計算量 $O(N ^ {2})$ という問題児な一面も兼ね備えています。

マージソート

分割統治の考え方に基づくアルゴリズムです。A[0:N)をソートすることを考えます。これをA[0:N/2)をソートする問題とA[N/2:N)をソートする問題に分割します。これを解いたら、この2つの数列をマージします。これは「A[0:N/2)とA[N/2:N)に対し、両者の先頭の値を比較して小さいほうを空の配列に加えて、元の配列から除去する」という操作を繰り返すことで可能です。言葉で聞くより図を見るほうがわかりやすいと思います。再帰の深さがおよそ $\log _ {2} N$ 、マージの計算量が各再帰深さについて合計 $O(N)$ なので、計算量は $O(N\log N)$ です。

このほか、挿入ソート： $O(N ^ {2})$ やヒープソート： $O(N\log N)$ などが有名だと思います。また、C++11以降のstd::sort()はイントロソートというアルゴリズムを使っています（規格として定まってはいませんが、そのように実装されていることが多いそうです。参考: sort - cpprefjp C++日本語リファレンス）。イントロソートは、クイックソートをベースに、再帰が深くなるとヒープソートに切り替えることによって最悪計算量が $O(N\log N)$ に改善されるというものです。PythonやJavaやRustではティムソートなんてのを使っているらしいです。

変なソートを作ろう

思い付きで変なソートを作ってみたいと思います。目標は「バブルソートよりはマシ」、つまり $O(N ^ {2})$ よりもよい計算量です。

平方分割バブルソート

上で述べた通り、バブルソートは計算量が $O(N ^ {2})$ と悪いのですが、これにマージソートの要素を加えて計算量を改善してみたいと思います。

Aをいくつかのブロックに分割することを考えます。何個に分割するかは全体を俯瞰してから適切なものを決めたいので、ここではB個のブロックに分割するものとします。すると、Aは長さN/Bの配列B個に分割されます。

次に各ブロックに対し、ブロック内の範囲でバブルソートします。この操作の計算量は全体で $O((N/B) ^ {2}\times B)=O(N ^ {2} / B)$ です。これによりソートされた長さN/Bの配列がB個できあがります。

最後にこれらの配列を1つの配列にマージします。B個の配列の先頭を走査し、最小の値をソート済み配列に追加し、元の配列から削除する、という一連の操作を繰り返します。配列の要素数の総和は当然Nですから、B要素の走査がN回行われることになり、この処理の計算量は $O(NB)$ です。

全体の計算量は $O(N ^{2}/B + NB)$ となりました。 $B$ はいくらにするのが最適でしょうか。大きくすると $NB$ が、小さくすると $N ^ {2}/B$ が支配的になってしまうので、その中間の $N ^ {2}/B=NB$ になるあたりで手を打つのがよさそうです。これを変形すると $B=\sqrt{N}$ になります。よって配列を長さ $\sqrt{N}$ の配列 $\sqrt{N}$ 個に分割して上の処理をすれば、ソート部分とマージ部分の計算量がどちらも $O(N\sqrt{N})$ になり、 $O(N ^ {2})$ から改善されます（もちろん通常の $O(N \log {N})$ のマージソートには劣ります）。

このように大きさ $\sqrt{N}$ のブロックに分割して計算量を落とす手法を平方分割と呼んだりします。競プロでもよく使う手法です。

例えばこんな問題が解けます

参考にしたのはこちらです。平方分割の練習をしようにも難しい問題ばかり、そんなお悩みに狙いを決めて手取り足取りのレクチャーです！ - CADDi Tech Blog

長さNの配列Aが与えられる。以下の2種の命令を順に処理せよ。

A[i]をxに変更せよ。
A[l:r)の最大値を答えよ。

上の命令は簡単で、即座にできます。下の命令は単純に計算すると $O(N)$ です。重たい処理に引っ張られるので、命令あたりの計算量は $O(N)$ です。

ここで配列を長さ $\sqrt{N}$ のブロック $\sqrt{N}$ 個に分割します。また事前に各ブロックごとに要素の最大値を計算し、メモしておきます。こうすると、範囲[l:r)はいくつかのブロックといくつかの端数に分割されます。ブロックの数は $\sqrt{N}$ 個以下、端数は右と左で各 $\sqrt{N}$ 未満ずつなので、それらの最大値を求めるための計算量は $O(\sqrt{N})$ に高速化されます（もちろん、変更命令の際にA[i]だけでなくA[i]が属するブロックの最大値も変更する必要があります。素直にブロック内を走査すればよく、 $O(\sqrt{N})$ で計算できます）。なお、この問題はsegment treeというデータ構造（配列の一点変更と区間最大値をどちらも $O(\log N)$ で処理できます、考え方はマージソートで使ったものに似ています）などを用いてもっと効率的に解くことができます。

ところで、このソートにおいて $\sqrt{N}$ 個の配列に分割した後、それらに対してこのソートを再帰的に適用すると、計算量は改善されるでしょうか。答えは否です。これは最終的に長さ $\sqrt{N}$ の配列 $\sqrt{N}$ 個をマージするパートの計算量が $O(N\sqrt{N})$ から落ちないためです。再帰的に処理する元気があるなら普通のマージソートを使ってください。

追記

先行研究（バブルソートか選択ソートかの違いでほぼ同じもの）が見つかったので紹介します、"ソート平方分割"などのワードで検索すると他にも出てくると思います：

平方分割でマージソートするやつ（O(N√N)なので実用性は無）https://t.co/NeZtktnHlB
— 🐟 (@0x3b800001) 2021年11月16日

最長増減部分列分割ソート

ソートの話をしすぎて、ソートをサボりたい気分になりました。配列からすでにソートされている部分を抜き出して、それらをマージするというのはどうでしょうか。例として

A={1,4,7,3,5,6,2,8}

を考えます。分割の際、先頭は必ず選ぶとして、それと同じか大きいものになっていくように選んでいきます。

A={1,4,7,3,5,6,2,8}

A={3,5,6,2}

A={2}

これで3個の配列{1,4,7,8} {3,5,6} {2}に分割されたのでこれを先ほどと同じ要領でマージします。

しかし、問題があります。こんな配列はどうでしょう。

A={8,7,6,5,4,3,2,1}

これは先頭が常に最大になるので8個の数列（{8}{7}{6}{5}{4}{3}{2}{1}）に分割されてしまいます。逆順ソートされた配列に対してはその要素数ぶんだけ長さ1の配列が作られるので、マージの計算量が $O(N^ {2})$ になってしまいます。処理が複雑なぶん、バブルソートにも負けてしまいそうです。

では、ソートされている部分を抜き出すのと逆順ソートされている部分を抜き出すのを交互にやったらどうでしょうか。こうすれば、先ほどの配列は{8}{7,6,5,4,3,2,1}に分割され、逆順ソートされているものをもとに戻してマージすればいい感じになります。最初の{1,4,7,3,5,6,2,8}もできる配列は3つのままです。

本当にこれでいいのでしょうか。答えは否です。例えば、

A={8,1,7,2,6,3,5,4}

は8つの配列に分割されてしまいます。（逆順ソートを先に抜き出すと若干改善しますが、ほぼ変わりません）

では、この考えかたでは $O(N ^ {2})$ から改善することはできないのでしょうか。もう少し頑張ってみます。

なぜ先ほどのやり方ではうまくいかないのでしょうか。それは、「目先のものに後先考えず飛びついているから」です。{8,1,7,2,6,3,5,4}がいい例で、先頭の8に飛びつかず、次の1が来るのを待って、7を我慢して、2を取って、という感じで、できれば{1,2,3,4}{8,7,6,5}に分割してもらいたいのです。ようは、できるだけ長い部分列を取ってきてほしいわけですが、これを実現するアルゴリズムがあります。

最長増加部分列

上の問題は最長増加部分列（Longest Increasing Subsequence、LIS）といい、競プロでも頻出です。これを解くための計算量は $O(N\log N)$ で、{8,1,7,2,6,3,5,4}から部分列{1,2,3,4}を抜き出すことができます。アルゴリズムの中身は動的計画法（Dynamic Programming）と呼ばれるものです。配列dpを用意して、Aを走査しながら

$dp [ i ]$ = 長さ $i$ の増加列の末尾の最小値

となるように更新していきます（今回はLISを配列から抜き出す必要があるので、末尾の最小値ではなく末尾の最小値が格納されているインデックス、としています）。

LISについてもう少し詳しく

例としてA={1,6,2,5,8,3,4,7}のLISを考えます。最初、dp={inf, inf, inf, inf, inf, inf, inf, inf}です（infは無限大の意）。また、今だけ1-indexにします（そのほうが理解が容易です）。

A[1]=1に注目します。これを使うと長さ1の増加列{1}ができます（増加列というのは違和感ですが）。dp[1]:=1と更新します。dp={1, inf, ..., inf}です。
A[2]=6です。先ほどの1と合わせて長さ2の増加列{1, 6}ができます。dp[2]:=6です。dp={1, 6, inf, ..., inf}です。
A[3]=2です。残念ながら長さ3の列{1,6,2}は増加列ではありません（6→2で減ってしまう）。ですが、最初の1と合わせて長さ2の増加列{1, 2}を作ることができます。これは{1, 6}よりも嬉しいです。なぜなら、6より大きい数よりも2より大きい数のほうが多いため、増加列を伸ばしやすくなるからです。dp[2]:=2、dp = {1, 2, inf, ..., inf} です。
A[4]=5です。さっきのdp[2]=2が効いて、長さ3の増加列{1, 2, 5}が作れます。dp[3]:=5、dp = {1, 2, 5, inf, ..., inf}です。
A[5]=8です。dp[4]:=8、dp = {1, 2, 5, 8, inf, ..., inf}です。
A[6]=3です。今度の更新はdp[3]:=3になります。5を3に書き換えることで後ろに続くことができる数の範囲を広げます。dp = {1, 2, 3, 8, inf, ..., inf}です。
A[7]=4です。dp[4]:=4とできます。dp = {1, 2, 3, 4, inf, ..., inf}です。
A[8]=7です。5. の終わりではdp={1, 2, 5, 8, ...}でしたがこれが6. と7. で{1, 2, 3, 4, ...}に変わったのが効いて、dp[5]:=7とできます。dp = {1, 2, 3, 4, 7, inf, ..., inf}でフィニッシュです。dpからinfでない部分を抜き出して、LISは{1, 2, 3, 4, 7}となります。

さて、dp[?]:=A[i]という更新をたくさんやってきましたが、この?の位置はどうやって決めればよいでしょうか。よく観察すると、この位置は「A[i]より大きい数が入っている場所のうち、最も左」です。これは二分探索で効率的に求められます。

（LISについてはこちらも合わせてご覧ください： LIS でも大活躍！ DP の配列使いまわしテクニックを特集 - Qiita）

二分探索とは

二分探索は単調増加/減少なものから何かを探索するときに使うアルゴリズムです。

今回は、「dp配列のうちA[i]より大きい数値が入っている場所の左端」を求めます。dp配列が以下のように単調増加になっているので二分探索が使えます。

$dp[1]\le dp[2]\le ... \le dp[N]$

まず、配列全体のど真ん中を見ます。もしその値がA[i]より大きいなら、更新すべき場所は配列のど真ん中それ自身か、それより左側のどこかです。逆に小さければ、更新すべき場所は配列のど真ん中より右側のはずです。これだけで考えるべき範囲が一気に半分になりました。次に、半分になった探索範囲のど真ん中を見て、それよりも右か左かを判定します。すると、全体の半分だった探索範囲をさらに半分にできます。これを繰り返すことで探索範囲を毎回半分に減らしていくことができ、いずれ1か所に特定できます。繰り返し回数はおよそ $\log _ {2} N$ 回で、Nに対してかなり小さい値になります。

このように、探索範囲を半分に狭めながら行う探索が二分探索です。

これを用いて配列の最長増加部分列と最長減少部分列のうち長いほうを順次抜き出していき、最後それらをマージする、これならいい感じに計算量が減ってくれそうです。

さて、最長増加部分列または最長減少部分列に分解していくと、結局いくつの配列ができあがるのでしょうか。これについて、興味深い定理があります。それがErdös Szekersの定理です。

a,b を正の整数とする。各項が相異なる長さ ab+1 の数列があるとき，以下の二つのうち少なくとも一つは必ず存在する。

長さ a+1 の部分列で，単調増加なもの

長さ b+1 の部分列で，単調減少なもの

(Erdös Szekersの定理とその証明 | 高校数学の美しい物語, https://manabitimes.jp/math/909)

ここでは証明には触れませんが（鳩の巣原理という競プロでもよく使う原理を使います。そこまで難しくないと思うのでぜひ引用元のページを読んでみてください）、要はどんな場合も現時点の配列の長さの平方根以上の長さの単調増加or減少部分列を取ることができると解釈できます。

この性質を知ったうえで、ある長さの配列が何個ぐらいの部分列に分割されるかを実際にプログラムで試してみます。具体的には漸化式

$f(0)=0,f(n) = f(n- \lfloor \sqrt{n} \rfloor ) +1$

を $n$ を変えて計算します。すると、

$f(1000)=63, f(10000)=199, f(10 ^ {6}) = 1999, f(10 ^ {10})=199999$

と、 $f(n)\approx 2\sqrt{n}$ になりそうです。

ざっくりとした理屈を考えてみましょう。まず $\lfloor\sqrt{n}\rfloor$ が気持ち悪いので $n=x ^ {2}$ として式を書き直すと、

$f(x ^ {2})=f(x ^ {2}-x)+1=f( (x-1/2) ^ {2}-1/4)+1\approx f( (x-1/2) ^ {2})+1$

となります。 $g(x)=f(x ^ {2})$ とおくと $g(x)=g(x-1/2)+1$ となり、これは一次関数 $g(x)=2x$ です。結果的に、

$f(n)=g(x)=2x=2\sqrt{n}=O(\sqrt{n} )$

になりそうです。

配列数が $O(\sqrt{N})$ なので、マージにかかる計算量は $O(N\sqrt{N})$ です。問題はLISを求めるパートの更新計算にかかる $\log {N}$ です。二分探索を行う回数が $O(N)$ とかになってくれないかなーと思ったのですが残念ながら $O(N\sqrt{N})$ になるようなので、こちらがボトルネックで $O(N\sqrt{N}\log{N})$ っぽいです。

二分探索回数の解析

配列がいくつに分割されるか考えた時とほぼ同様の手法です。その時点での配列長 $n$ における二分探索の実行回数 $S(n)$ は、

$S(0)=0,S(n)=n+S(n-\lfloor\sqrt{n}\rfloor)$

であると考えられます。 $n=x ^ {2}$ として書き直すと

$S(x ^ {2})\approx x ^ {2}+S( (x-1/2) ^ {2}) =x ^ {2}+(x-1/2) ^ {2}+S( (x-1) ^ {2})=...=\sum _ {i=1} ^ {2x} (i/2) ^ {2}=(\sum _ {i=1} ^ {2x} i ^ {2}) /4$

で、自然数の二乗の和の公式

$\sum _ {i=1} ^ {x} i ^ {2}=x(x-1)(2x-1)/6$

より、

$S(n)\approx (16x ^ {3}/6) / 4=2x ^ {3}/3=O(x ^ {3})=O(n\sqrt{n})$ です。実際にプログラムで検証してみると、これに近い結果が得られます。

実験

手元のノートパソコンで、要素数N={10000, 50000, 100000, 500000}、0以上 2³⁰ 未満のランダムな値の配列各100個に対して各種ソートを実行し、実行時間を計測してみます。せっかく計算量を評価したので、変なソートたちがバブルソートより速いこと（そして標準のソートより数段遅いこと）、あとできれば変なソート達の計算量の差が確認できると嬉しいです。

というわけでプログラムを書きました。200行近いソースコードをここに貼ると大変なことになるのでもう少し下にリンクを貼っています。

下表に平均処理時間[ms]を示します。xと書いてある部分のデータは膨大な時間がかかるのでとっていません。処理時間がきわめて短いものは範囲で示しています。

N	標準ソート	バブルソート	平方分割	部分列分解
1e4	0~1	118	1~4	12
5e4	1~4	3915	27	158
1e5	4~8	x	74	463
5e5	34	x	740	5732

この表から分かる通り、C++標準のソートは速いしバブルソートは使い物になりません。平方分割によってバブルソートは格段に速くなりますが標準のソートよりだいぶ遅いです。部分列分解は平方分割よりさらに数倍遅く、 $\log$ や実装の重さが響いている感じでしょうか。

ソースコード

リンクはこちら（N=100や1000のデータもとるようになっているのですが、軽すぎて無意味でした）。自然な実装にしてあるので、「おいおい何やこの無駄な処理は」と思う方は我慢してください。Intel Core i7-7500U(3.5GHz)、gcc8.1.0、C++17、最適化オプションO2です。

strangesort.cpp - Google ドライブ

おわりに

使い物にならないソートを作り、その性能を確かめました。皆さんも面白いソートを作ってみてはいかがでしょうか。

点Nの軌跡

競プロの話または日記