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

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 よりもだいたい実装が軽くなるのでオススメです．

a, b を互いに素な整数としたとき，方程式 a x + b y = 1 の整数解を求める手続きは拡張ユークリッド互除法（extended Euclidean algorithm）として知られています．今回はこれに関連する話題を紹介します．

非再帰のユークリッド互除法

http://topcoder.g.hatena.ne.jp/iwiwi/20130105/1357363348 で Mi_Sawa さんが紹介している非再帰版のユークリッド互除法（以下の手続き）を導出するところからはじめます．

LL extGCD(LL a, LL b, LL &x, LL &y) {
  for (LL u = y = 1, v = x = 0; a; ) {
    LL q = b / a;
    swap(x -= q * u, u);
    swap(y -= q * v, v);
    swap(b -= q * a, a);
  }
  return b;
}

まず次の行列を考えます．

A' = |a 1 0|
     |b 0 1|

（プライムは後の都合による）．すると，x = [-1, a, b] に対して A' x = 0 が成立します．そこで，この等式を崩さないように A' を行基本変形していきます．整数行列に対する行基本変形は

• ある行に別の行の整数倍を加える
• ある行の符号をかえる
• ある行と別の行を入れ替える

という操作の繰り返しでした．これらは全て行列の左から 2×2 行列 S をかけることで実現できます．A' x = 0 の左から S をかけると S A' x = 0 なので，行基本変形は「x を解とする」ことを保ちます．また，基本変形がすべて可逆なので S は正則であり，S A' x = 0 から A' x = 0 が導けることも注意しておきます．

さて，行基本変形の目的地として，一番左の列が簡単な形になるようにしましょう．基本変形の操作でユークリッドの互除法がシミュレートできることに注意すると，基本変形で

S A' = |g x y|
       |0 u v|

までいけることがわかります．ここで g = gcd(a,b) です．これが [-1,a,b] を解にもつ，ということを書くと

  a x + b y = g
  a u + b v = 0

となって不定方程式の解がもとまっていることがわかります．以上の手続きをコードにすると冒頭のものになり，非再帰の拡張互除法が導出できました．

...

Blankinship's algorithm

上の形でのユークリッドの互除法の理解は [Blankinship 1963] によります．一般に，不定方程式

　a_1 x_1 + ... + a_n x_n = gcd(a_1, ..., a_n)

を解くためには行列 A' = [A, I] （A は a が縦に並んだ行列，右は辻褄の合う単位行列）について基本変形を行えばよい，ということを Blankinship は指摘しました．この意味で冒頭の非再帰版の拡張互除法は Blankinship's algorithm とも呼ばれます．

なお，この手続きによって実は不定方程式の一般解も求まっている，ということが [Morito-Salkin 1979] によって指摘されています．上の２変数の例でいえば，すべての解は [x,y] + k [u,v] で書ける，ということですね．

typedef long long LL;
LL extGCD(LL a[], LL x[], int n) {
  LL A[n][n+1];
  memset(A, 0, sizeof(A));
  for (int i = 0; i < n; ++i) {
    A[i][n] = a[i];
    A[i][i] = 1;
  }
  while (1) {
    int k = -1; // k = nonzero argmin A[*][n]
    for (int i = 0; i < n; ++i) 
      if (A[i][n]) 
        if (k == -1 || abs(A[k][n]) > abs(A[i][n])) k = i;
    bool fin = true;
    for (int i = 0; i < n; ++i) {
      if (i == k || A[i][n] == 0) continue;
      fin = false;
      LL q = A[i][n] / A[k][n];
      for (int j = 0; j <= n; ++j) 
        A[i][j] -= q * A[k][j];
    }
    if (fin) {
      for (int j = 0; j < n; ++j)  // A[k] = x; a x = g (particular solution)
        x[j] = A[k][j];           // A[!k] = u; a u = 0 (homegeneous solution)
      return A[k][n];
    }
  }
}
int main() {
  int n = 3;
  LL x[n], a[n] = {12,15,10};
  LL g = extGCD(a, x, n);
  for (int i = 0; i < n; ++i)
    cout << x[i] << " "; // -2 1 1
}

線形方程式の求解として見る

※コードここまで

上の手続きの別の見方として，線型方程式 [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): 線形代数汎論, 朝倉書店.

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

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

劣微分

まずは連続の場合から復習します．関数 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):

今回は乗算型重み更新（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.