週刊 spaghetti_source

できるだけ毎週土曜日更新．次回更新は2月16日．
ネタ切れ深刻なので記事リクエスト募集中．twitterないしはコメント欄まで．
貼っているコードはあんまり検証していないので，信頼性はアレです．アレ．
twitter: @tmaehara フォローはお気軽に．

2013-02-09

Rolling Hashing

18:57

文字列系の問題に対するインチキ・データ構造を紹介します．

Rolling Hashing

文字列 s = s[0,n) に対してハッシュ値を次のように対応させます：

　h(s) := (s[n-1] + p s[n-2] + ... + p^{n-1} s[0]) mod M

つまり文字列を p 進数表示だと思って mod M で値をとったものです．Rolling Hashing とは，すべてのプレフィックス（もしくはサフィックス）に対してこの値を計算することをいいます：

　phash[k] := h( s[0,k) ) for k = 0, ..., n

ここで p と M は互いに素な数を選びます．以下では p = 10^9+7，M = 2^64 に選びます．

Rolling Hashingで得られるテーブルを Rolling Hash Table ということがあります．また，Prefix Hash とか（これは同名があってヤヤコシイ），Sequence Hashと呼ばれることもあります．

部分文字列に対するハッシュ値計算

すべてのプレフィクスに対してハッシュ値を計算するのは普通に前から計算するだけです．そのテーブルがあれば，任意の部分文字列に対するハッシュも h(s[i,j)) = h(s[0,j)) - h(s[0,i)) * p^{j-i} として計算できます．p^{k} をすべての k について前計算しておけば h(s[i,j)) も O(1) です．

  static const ULL p = 1000000007ULL;
  string s;
  int n;
  vector<ULL> pow, phash;
  RollingHash(string s) : s(s), n(s.size()), pow(n+1), phash(n+1) {
    pow[0] = 1; phash[0] = 0;
    REP(i, n) {
      phash[i+1] = s[i] + phash[i] * p;
      pow[i+1] = pow[i] * p;
    }
  }

...

で，何ができるの？

ハッシュ値の衝突を無視すると，部分文字列同士の一致比較が O(1) できます．

いわゆるsuffix arrayが部分文字列の大小比較が O(1) でできる構造なので，rolling hash tableは劣化版suffix arrayだと思って良いかと思います．出来ることも基本的にニア・イコール．メリットは実装の簡素さにあります．

文字列検索（Rabin-Karp）

パターンの hash を計算して，テキストの前から順番に hash 値を比較するだけ（Rolling Hash で前計算しておく）．この手続きはRabin-Karpとして文字列検索では良く知られています．普通の文字列検索ではKnuth-Morris-PrattやBoyer-Mooreがあるのであまり使いませんが，基本事項です．

  ULL hash(string t) {
    ULL h = 0;
    REP(i, t.size()) h = t[i] + h * p;
    return h;
  }
  int find(string t) {
    int w = t.size(), count = 0;
    ULL h = hash(t);
    REP(i, n-w+1) if (hash(i, i+w) == h) {
      ++count; /* found */ 
    }
    return count;
  }

最長共通接頭語：lcp(i,j) の計算

LCPの計算は，普通なら suffix array を作る際に O(n) で追加のテーブルを作っておいて RMQ を使って O(1) で計算するんだ！というところですが，rolling hashing を使っても O(log n) で計算できます．すなわち，lcp(i,j) は s[i,i+k) == s[j,j+k) を満たす最大の k なので，これを二分探索します．

  int lcp(int i, int j) {
    int l = 0, r = min(n-i, n-j)+1;
    while (l + 1 < r) {
      int m = (l + r) / 2;
      if (hash(i,i+m) == hash(j,j+m)) l = m;
      else                            r = m;
    }
    return l;
  }

Suffix Arrayの計算

k = lcp(i,j) とすると，s[i,n) と s[j,n) の大小関係は s[i+k] と s[j+k] の大小関係に一致します．なので，コレを operator < に採用して sort するだけで suffix array が作れます．計算量は O(n (log n)^2)．定数が小さいのがポイント．

struct RollingHash {
  // ...
  bool operator()(int i, int j) { 
    int k = lcp(i, j);
    return i+k >= n ? true : j+k >= n ? false : s[i+k] <= s[j+k];
  }
  // ...
};

  RollingHash R(s);
  int n = s.size();
  vector<int> sa(n);
  REP(i,n) sa[i] = i;
  sort(ALL(sa), R);

最長回文

部分文字列 t のうち t == rev(t) が成り立つ最長のものを求める問題を考えます．

簡単のため奇数長の回文に限ります．すると，長さ k の回文があれば長さ k-2 の回文もあるので，「最長回文の長さ」を二分探索できるようになります（偶数長への対応は間にダミーを挟むなどで可能）．

長さ k の回文があるか，という質問は，O(n) で長さ k の部分列のハッシュを全部表に登録しておいて，reverse(s) の長さ k の部分列を O(n) で見てハッシュで一致するものがあるかを調べます．全体で O(n log n)．

最長k出現部分列

部分文字列の中で k 回以上出現するもののうち，最も長いものを求める問題を考えます．

問題を反転して，長さ m の文字列のなかで最も出現回数が多いものの回数を f(m) とおくと，f(m) は m に関する単調減少関数なので，m に関する二分探索で f(m) >= k なる最大の m を求めるのが全体的方針．

f(m) の計算は長さ m の部分文字列のハッシュ値を全部テーブルに突っ込みつつ頻度の max をとるだけ．

...

このように，割と単純な手続きで色々と高効率なものが作れます．特に suffix array よりもだいたい実装が軽くなるのでオススメです．

コメントを書く

上の手続きの別の見方として，線型方程式 [x,y] A = c を解く問題だと考えます．右辺は拡張互除法では g でしたが，一般的に c にしておきます．普通の線型代数で方程式 x A = c を解こうと思ったら，A を適当に基本変形してくことになるので，それに倣って基本変形していきます．ただし，今 A はサイズ的に左からの積しか意味が無いので，左から操作します．左からの積は行基本変形に対応していたことに注意すると

A = |a|
    |b|

に適当に左から行列 S をかける，ということなので

S A = |g|
      |0|

となって，x S^{-1} = [s,t] とおけば結局 x A = s g = c となって，c が g で割り切れるときのみ整数解をもち，「s = 1/g, t = 任意」がその解全てであることがわかります．[s,t] から x を求めるには [x,y] = S [s,t] なので S が求まれば十分です．これは学部の線形代数の講義を思い出すと，A を A' = [A I] と拡張しておいて A 側に基本変形を行なって I に伝播させればよい，ということでした．この手続きは Blankinship の手続きと一致しています．

※学部の線型代数でやることは，行列 A の逆行列を求めるために [A I] を考えて基本変形で [I S] へと変形すると S = A^{-1} になる，というものだったはずです．今は A の逆行列まではいけないので「S A = 簡単な形」とする S をあらわに求める手続きになっています．

単因子標準形

さらに拡張します．A を整数係数の行列，b を整数ベクトルとして，線型方程式 x A = b を満たす x を求める問題を考えます．今度は A は行基本変形も列基本変形も意味があるので両方やります．どこまで A がいけるか，ということに答えるものが単因子標準形です．

定理． 整数を成分とする行列 A にたいし，あるユニモジュラ行列 S, T が存在して

S A T = diag(g[1], ..., g[r])　（サイズのあわないところはゼロ）
g[1] = gcd(A[*]), g[i] は g[i+1] を割り切る

例で見ましょう．

 |2 6| → |2  6| → |2 2| → |2 0|
 |4 8|    |0 -4|    |0 4|    |0 4|

つまり，上の A = [2,6;4,8] という行列は S A T = diag(2,4) となる，ということです．方程式 x A = b を解きたければ A を単因子標準系にして x S^{-1} = y とおいて y (S A T) = b T とし，y(i) = (B T)(i) / g(i) が解なのでこれを x = S y で戻せば良い，ということになります．これで一般の整数係数 A x = b が解けるようになりました．

実装では，基本変形は絶対値最小の要素を探してそれで行と列を掃き出すこと有限回で終わります（max |a*| が単調減少する）．これは１列の場合は Blankinship の手続きと一致するので，結局 Blankinship の拡張互除法は単因子標準形を求める掃き出し法の特殊版であったと理解できます．

一般に，ここまで述べたことはすべてユークリッド整域（ユークリッド互除法が定義できる領域）で成り立ちます．ユークリッド整域の例には他にも Z，Z[√-1], Z[√-2], Z[ω], K[x] などがあります．複素整数環 Z[√-1] や多項式環 K[x] の gcd は一度は書いておきたいかもしれません．

参考文献

W. A. Blankinship (1963): A new version of the Euclidean algorithm, American Math. Monthly, no.7, pp.742-745.
S. Morito, and H. M. Salkin (1979), Finding the general solution of a linear diophantine equation, Acta Informatica, vol.13, pp.379-382.
伊理正夫 (2009): 線形代数汎論, 朝倉書店.

コメントを書く

2013-01-12

subgradient representation of matroid intersection theorem

21:03

今回は劣微分によるマトロイド交差定理の表現を紹介します．この形でマトロイド交差定理を説明している文献を見たことがありませんが，従来表現から簡単に出るので特に新しくはありません．ただ，自分の肌には良くあう理解です．

以下では，まず劣微分について簡単に復習し，その後，マトロイド交差定理を表現します．

劣微分

まずは連続の場合から復習します．関数 f: V→R の点xにおける劣微分とは，点xでfを下から支持する超平面のことです．式で書くと次のようになります：

　∂f := { p \in V^* : (p,y-x) ≦f(y)-f(x) for all y }

(・,・)はペアリング（内積）です．任意の点で劣微分がとれることと，凸であることはほとんど同値（端点だけ例外）となります．病的な境界の関数を考えると違いが出てくるのですが，実用上は「凸解析＝劣微分がとれる関数クラスの理論」と思ってほぼ正しいです．

劣微分と最適化の繋がりで最も重要なのは次の命題です：

命題． 0∈∂f(x) であることと，x が f の大域最適解であることは同値．

（∵ z軸に並行な接平面がとれる＝自分より下が居ない．//）

このことは「0∈∂f(x)」が「x が最適である証拠」だということを主張しています．逆に，最適の証拠がとれない場合には，劣微分を用いて最適解方向へ x を更新することが可能です：

命題 [see: Boyd 2003]． 任意の 0≠p∈∂f(x) について，あるh > 0が存在して，d(x*, x - h p) < d(x*, x)．ここで x* はある大域最適解．

（∵ d(x*, x-hp) を計算する．やるだけなので省略．//）

※微分可能の場合によく使う勾配方向への更新では必ず目的関数値が減少しますが，劣微分方向への更新では減少しない例が作れます．それでも，x 自体は最適解へと近づきます．これは「劣勾配法の強収束性」と呼ばれています．

これらのことから，標語的に

　凸関数 <=> 劣微分存在 <=> 大域最適が計算可

という関係が成り立ちます．

さて，劣微分に関する重要な性質として，次の定理（劣微分の和公式）があります．

定理．[Moreau-Rockafeller; see: Winfried 2007] f, g を凸関数とする．弱い仮定のもとで ∂(f + g) = ∂f + ∂g．

一般の f, g で ⊇ 側の包含関係は従いますが，逆側が従うのが凸関数の良いところです．証明は凸集合の分離定理（Hahn-Banach）を用います．

※弱い仮定は「エピグラフに対して凸分離定理が成り立つ」程度の仮定で，通常は成立すると思ってよい．

マトロイド

さて，ここからが本題です．マトロイドに関して劣微分の議論を導入します．関数を考えるので，マトロイドよりもM凸関数を考えたほうが自然です．関数 f がM凸であるとは，任意の x, y と i ∈ supp(x-y) について，ある j ∈ supp(y-x) が存在して

　f(x) + f(y) ≧ f(x-ei+ej) + f(y+ei-ej)

が成り立つことをいいます．たとえばマトロイドに対して f を次のように定義するとM凸関数になります．

f(x) = -|x| if x 独立, +∞ if x 従属

※一般に，M凸関数の実行定義域（値が有限をとる点）全体はマトロイド（の拡張）となります．マトロイドはM凸関数の実行定義域のことだ，と思っておいて大丈夫です．

M凸関数 f の劣微分はどうなるでしょうか．結論から言うと，これは

　∂f(x) = { p ∈ R^n | p(j) - p(i) ≦ f(x-i+j) - f(x) }

となります．最適性条件 0 ∈ ∂f(x) を書き下すと，上の式で p = 0 とおいて

　0 ≦ f(x-i+j) - f(x) (for all i, j)

となります．これはM凸関数についてよく知られている最適性条件です（実際は：最適性条件から劣微分を出すほうが易しい）．

最適でない場合の更新ですが，連続関数と違い，劣微分方向に進むことはできません．そのかわり，劣微分方向で整数格子上の一番近い点へと丸めることで，いわゆる「マトロイド上の貪欲法」が出てきます（減り幅最大の方向に進む）．つまりマトロイド貪欲は劣勾配法の一種と言えます．

続いて，マトロイド交差問題は次のように表現できます：f, g M凸関数について

　min f(x) + g(x)

x が f, g どちらかの実行定義域（独立集合）を外れると全体が無限大に飛ぶので，共通実行定義域（＝共通独立集合）に入ることが x の条件です．このことから，確かにマトロイド交差問題を表していることがわかります．

この問題の最適性条件を見るために ∂(f + g) について考えます．実は，マトロイド交差定理は，Moreau-Rockafellerの定理の離散版と言えます：

定理 [Matroid Intersection; see: Murota 2003]． f, g をM凸関数とする．∂(f + g) = ∂f + ∂g．

（∵Frank's weight splittingを使う）

※M凸でない関数和では成り立ちません．例えば

　f(x,y) = |x + y - 1| (… M凸)

　g(x,y) = |x - y| (… L凸)

はどちらも連続の意味では凸関数ですが，離散に制限すると和公式が成立しません．実際 f + g は (0,0) で最小値 1 を達成するので 0 が(f + g)の劣微分に入るはずですが，f の (0,0) での劣微分は (-1,-1) の一点のみ，g の (0,0) での劣微分は { (t,-t) | -1 ≦ t ≦ 1 } なので，和でゼロは作れません．

...

劣微分からアルゴリズムへ

和公式が成り立つことを信じると，マトロイド交差アルゴリズムが自然に導出できます．最適性条件 0 ∈ ∂(f + g) = ∂f + ∂g は，ある p が存在して p ∈ ∂f, -p ∈∂g となることと同値です．これをM凸の劣微分の式とあわせると

　p(j) - p(i) ≦ f(x-i+j) - f(x)

　p(i) - p(j) ≦ g(x-i+j) - g(x)

となります．ポイントは，この不等式系は最短路問題の双対 [Cormen-Leiserson-Rivest-Stein 2009] となっていることです．すなわち，グラフ G を「辺 (i,j) のコストが f(x-i+j) - f(x)，(j,i) のコストが g(x-i+j) - g(x) となるグラフ（＝交換グラフ）」とすると，p が存在することはグラフが負サイクルをもたないことと同値となり，負サイクルがないときは p として最短路長がとれます（適当に p(s) = 0としてよい）．

逆に，負サイクルがあった場合，そのサイクルから劣微分の元を構成することができます．すなわち，負サイクルキャンセルアルゴリズムは，適当な劣微分（負サイクル）による劣勾配法とみなせます．

...

まとめ

マトロイド交差定理は ∂(f + g) = ∂f + ∂g と表現可能．
※連続凸なら常に成り立つが離散凸では成り立つとは限らない．
f+g の最小値を求めるために ∂(f+g) = ∂f+∂g = 0 とすると，交差問題の最適性条件（交換グラフに負閉路なし）が出る．
負閉路がある場合にそれにそってキャンセルするのは，劣勾配法の一種．

参考文献

S. Boyd (2003): "Subgradient methods", Notes for EE392o, Stanford University, Autumn, 2003.
S. Winfried (2007): "Nonsmooth Analysis", Springer Berlin Heidelberg.
K. Murota (2003): "Discrete Convex Analysis", Society for Industrial and Applied Matheamtics.
T. H. Cormen, C. E. Leiserson, R. L. Rivest, and C. Stein (2009): "Introduction to Algorithm", The MIT Press, 3rd ed.

S. Boyd (2003):

コメントを書く

2013-01-05

multiplicative weights update method

15:59

今回は乗算型重み更新（multiplicative weights update method）について説明します．いろんな問題に適用できる手法ではありますが，基本的には「連続最適化問題の近似解を得る手法」なので，離散最適化問題・厳密解を重視する競技プログラミング的にはあまり役に立たなそうなネタです．ただし，実装は超容易なので，覚えておくと得することはあるかもしれません．

...

乗算型重み更新は双対勾配法の一種 [Baes-Burgisser 2010] で，様々な分野で何度も再発明されているテクニックです．特に機械学習分野では標準的な技法だったようですが，最適化分野で注目が高まったのは [Arora-Hazan-Kale 2007] による「sparsest cutのSDP緩和を高速に解く」という話からでした．この結果以降，サーベイ論文 [Arora-Hazan-Kale 2012] （2010年にプレプリントとして公開）を参考にした Lecture Note が多くの最適化分野の人によって書かれており，現在では標準的手法の１つになった感があります（see: Google検索： multiplicative weights, lecture note）．

なお，より詳しい歴史的背景は [Arota-Hazan-Kale 2012]を参照してください．今回の解説もこのサーベイ論文を参考にしています．

さて，以下では乗算型重み更新の説明を行いますが，これは「幅広い問題に対するアルゴリズムの設計方針」みたいなものなので，抽象化して説明するのが大変です．自分が思うに，乗算型重み更新を理解するには，いろんな例題に対して適用してみるのが良さそうなので，今回は典型的な２つの問題に対して適用してみます．照明のパターンが重要なので，興味のある方は是非紙とペンを持って計算をフォローしてみてください．

...

問題１．オンライン意思決定問題

次の状況設定を考えましょう．

n 個の判別機がある．
質問が繰り返し与えられる．各質問に対して
- 判別機は独立に答え（Yes / No）を出力する．
- 自分も何かしらの答えを出力する．
- 質問の正答が教えられる．

この設定で，トータルの自分の正答数をできるだけ多くするのが目的となります．どのように各判別機の結果を参考にするかが戦略（意思決定）にあたります．乗算型重み更新では，この問題に対して次のアルゴリズムを設計します．

重み p := (1, ..., 1) ∈ R+^n と初期設定する．
各質問に対して以下を実行する：
- 重みに比例する確率で各判別機の回答を採用し，出力する．
- 間違えた判別機の重みを exp(-h) 倍する（h は小さな正の数）．

イメージとしては「たくさん間違えるとどんどん重みが軽くなり，答えが採用されなくなる」という感じです．重みの変化のさせ方には様々なバリエーションがあると思いますが，上のように exp 倍で重みを変化させるテクニックのことを乗算型重み更新といいます．

※exp(-h) 倍のかわりに (1 - h) 倍するものも乗算型重み更新と呼びますが，exp の展開を考えると本質的に同値なので区別しません．一応 exp(-h) 型更新を「Hedge algorithm」といい，(1 - h) 型更新を「multiplicative weights update」と呼んで区別することもあるようです．

乗算型重み更新で重要なことは，結果の評価にあります．今回の場合次の結果を証明することができます：

定理． 自分が間違える回数は，判別機の中の最多間違い回数 + log n / h + O(h) 倍以下．

証明． k 番目の質問のときの重みを p(k) とし，k 番目の質問に対する各判別機の成果を u(k) とします（正答なら u(k)_j = 0，誤答なら u(k)_j = 1）．判別機 j の k 回目までの間違え回数が Σ[l=1:k] u(l)_j となることに注意してください．

ポテンシャル関数 F(k) := p(k)_1 + ... + p(k)_n とおいて，質問に対してどのように変化するかを評価します．

　F(k+1) = Σp(k)_j exp(-h u(k)_j)

ですが，右辺を h が十分小さいとして展開すると

　F(k+1) = Σp(k)_j ( 1 - h u(k)_j ) = F(k) - h (p(k), u(k))

となります．ここで (p(k), u(k)) / F(k) は k 番目の質問に対して自分が間違える確率に等しいので，それを m(k) という記号であらわすことにすると

　F(k+1) = F(k) - h F(k) m(k) = F(k) ( 1 - h m(k) ) = F(k) exp(-h m(k))

となります．ただし最後の等号では再び exp の近似展開を使っています．この式を繰り返し代入すると

　F(k+1) = F(1) exp(-h Σm(l)) = exp(log n - h Σm(l))

となりますが，一方で，F の形から，任意の j について

　F(k+1) ≧ p(k+1)_j = exp(-h Σ[l=1:k] u(l)_j )

という評価も成立します．これらをあわせて

　Σm(l) ≦ log n / h + Σ[l=1:k] u(l)_j

となります．誤差のオーダを見積もると + O(h) であることもわかります．//

...

証明のポイントは２点ありました．

F(k) と F(k+1) の差の評価で「自分の損失（間違え確率）」が出ること．
p(k)_j の評価で「各判別機の損失」が出ること．

特に各評価が exp の肩に乗っていることから，全体としてO(log n)の評価できることに強みがあります．全体としてO(n)の誤差でいいなら普通の -h 的更新でも可能です．

※上のように確率的にせず，決定的に行なっても類似の評価ができますが，評価の質が劣化します．これは「完全ランダム戦略」で1/2正答率が保証されることと関係があります．

---

問題２．二人ゼロ和ゲームの鞍点

行列 A に対し，min[x] max[p] p A x を満たす確率ベクトル x, p を求める問題を二人ゼロ和ゲームの鞍点問題といいます．ゲーム理論で有名なNash均衡の存在定理は min[x] max[p] = max[p] min[x] が成り立つ，というものでした．

この問題は，普通は線型計画問題に帰着して解きますが，乗算型重み更新でも解くことができます．

p = (1, ..., 1) と初期化する．
c := p A について min[x] c x を満たす確率ベクトル x を求める．
p_j = exp( h (A x)_j ) と更新する．

定理． p(k)/Σp(k)_j := q(k) とおく．このとき Σq(l) A x(l) / T は鞍点の値の log n / h T 近似である．

証明． 先程と同じく F(k) := p(k)_1 + ... + p(k)_m と置いて評価してみます．まず F(k) の漸化式は

　F(k+1) = F(k) + h (p(k) A x(k)) = F(k) (1 + h q(k) A x(k)) = F(k) exp(h q(k) A x(k))

となります．漸化式を繰り返し適用して

　F(k+1) ≦ exp(log n + h Σq(k) A x(k))

となりますが，一方で

　F(k+1) ≧ p(k+1) = exp(h Σ(A x(l))_j)

という評価も得られるので，x(l) の平均を単に x とおいて整理すると

　(A x)_j ≦ log n / h T + 1/T Σq(l) A x(l)

となります．これが任意の j について成り立つので，任意の確率ベクトルと内積して

　q A x - log n / h T + ≦ 1/T Σq(l) A x(l)

と整理しておきます．いま，左辺を q で最大化すると max[q] q A x ≧ min[x] max[q] q A x = OPT であり，右辺は x のとり方から q(k) A x(k) = max[x] q(k) A x ≦ min[q] max[x] q A x = OPT となるので

　OPT - log n / h T + ≦ 1/T Σq(l) A x(l) ≦ OPT

となって，1/T Σq(k) A x(k) が鞍点の log n / h T 近似になっていることが示されました．//

ここからもう少し頑張ると Σx(l)/T と Σp(l)/T が鞍点の近似解になっていることがわかります．

実装例

２つ目の例題：二人ゼロ和ゲームの鞍点，の実装を示します．実装が恐ろしくシンプルなことに注目してください．

def mult_weight():
  epsilon = 0.1;
  x = [1.0 for i in range(n)];
  y = [0.0 for j in range(m)];
  T = int(n / (epsilon*epsilon) );
  for t in range(T):
    c = [sum(x[i] * a[i][j] for i in range(n)) for j in range(m)];
    j = max(zip(c, itertools.count()))[1];
    y[j] += 1.0;
    x = [x[i] * (1 - epsilon * a[i][j]) for i in range(n)];
  return (normalize(x), normalize(y));
def normalize(x):
  return map(lambda t: t / sum(x), x);

参考文献

S. Arora, E. Hazan, and S. Kale (2007): A combinatorial, primal-dual approach to semidefinite programs, Proceedings of the 39 annual ACM symposium on Theory on computing, pp.227-236.
M. Baes and M. Burgisser (2010): Hedge algorithm and subgradient methods, Optimization Online.
S. Arora, E. Hazan, and S. Kale (2012): Multiplicative weights update method: a meta-algorithm and its applications (Reserarch Survey), Theory of Computing Journal, vol.8, no.6, pp.121-164.

コメントを書く

2012-12-16

Trie for Fixed Universe Set

19:14

要素の追加・削除・存在判定に答えるためのデータ構造として平衡二分探索木がありますが，与えられるものが有限集合：{0, ..., 2^M-1 (=U)} の元である，という情報がある場合はより高速なデータ構造が設計できます．この手の問題は Fixed Universe Problem などと呼ばれています．

この手の問題の最も基本的な方針は，数字の上位ビットと下位ビットに分けて，上位ビットをキーとする再帰構造に下位ビットを突っ込む，というものです．van Emde Boas木として知られている構造はMビットユニバースに対する上位ビットサイズを H(M) = M/2 に選んだデータ構造で，木の深さが log(M) = log log(U) となります．ところが，vEB木はその計算量インパクトと比べてほとんど実用はされていません．これはvEBの計算量ゲイン：log(U) vs loglog(U) が効いてくるには十分大きなユニバースでないといけないのに対し，そのような大きなデータの処理では，データアクセスのオーバーヘッド（ファイルの読み出し等）が支配的になってしまい，もっと単純でデータ配置に即したB木などの構造のほうがトータルで早くなるからだと思っています（see: 参考文献）．

※vEB拡張系（＝log log U系）で実測すると結構早い，というのもどこかで聞いた気はしますが，忘れてしまいました．実用しているという話は聞かないので，多分そういうことだと思います．

...

ここでは最もシンプルな fixed universe problem の実装として，上位ビットを固定で H ビットとする構造の実装を示します．これは数を H ビットごとで切った文字列と見たときの Trie です．

切るビット数 H の設定ですが，子が全部で 2^H 子あり，どの子が居るか居ないかの状態が 2^2^H 状態あるので，2^2^H = 2^31 から H = 7 程度にします．このようにすると，たとえば「自分の子の中で空でない最小の子」とかがビット演算によって非常に高速に求まります（＝最小値探索・successor探索なども高速）．このときの木の深さは 31/7 < 5 より，深さ高々 5 が従います．これは 2^31 に対するvEB木と同じ深さ： log(32) = 5 なので，この範囲では完全に有利になっています．

実測

普通の平衡二分探索木（std::set）と insert/remove/find の比較を行います．1,000,000個の乱数をランダムに挿入・削除・検索する，という試行を10回行ったときの合計時間を見ています（※削除・検索が自明にならないように，乱数範囲は適当に設定）．この結果から，insertは19倍，removeは4倍，findは9倍ほど早くなっていることがわかります．

name	insert	remove	find
trie	0.129	0.099	0.136
std::set	2.426	0.398	1.279

実装例

const int B = 7;
const int MAXBIT = 32;
struct IntegralTree {
  struct Node {
    Node *child[(1<<B)];
    int S;
    Node() : S(0) { fill(child, child+(1<<B), (Node*)0); }
  } root;
  void insert(Node *n, int x, int R) {
    int h = (x >> (R-B)) & ((1 << B)-1);
    n->S |= (1 << h);
    if (R > B) {
      if (!n->child[h]) n->child[h] = new Node;
      insert(n->child[h], x, R-B);
    }
  }
  void remove(Node *n, int x, int R) {
    int h = (x >> (R-B)) & ((1 << B)-1);
    if (!(n->S & (1 << h))) return;
    bool cond = true;
    if (R > B) {
      remove(n->child[h], x, R-B);
      if (n->child[h]->S != 0) cond = false;
    }
    if (cond) n->S &= ~(1 << h);
  }
  bool find(Node *n, int x, int R) {
    int h = (x >> (R-B)) & ((1 << B)-1);
    if (!(n->S & (1 << h))) return false;
    if (R > B) return find(n->child[h], x, R-B);
    else       return true;
  }
  void insert(int x) { insert(&root, x, MAXBIT); }
  void remove(int x) { remove(&root, x, MAXBIT); }
  bool find(int x) { return find(&root, x, MAXBIT); } 
};

参考文献

T. H. Cormen, C. E. Leiserson, R. L. Rivest, and C. Stein (2009): "Introduction to Algorithms", The MIT Press; 3rd ed.
P. Patel and D. Garg (2012): Comparison of Advance Tree Data Structures, International Journal on Computer Applications, vol.41, no.2.
J. Levandoski, D. Lomet, and S. Sengupta (2013): "The Bw-Tree: A B-tree for New Hardware", In: 2013 IEEE 29th International Conference on Data Engineering (ICDE).

コメントを書く

週刊 spaghetti_source

2013-02-09

Rolling Hashing

2013-01-26

Euclidean algorithm (Blankinship algorithm)

2013-01-12

subgradient representation of matroid intersection theorem

2013-01-05

multiplicative weights update method

2012-12-16

Trie for Fixed Universe Set