二重確率行列を求める - らんらん技術日記

（ Bilateral Solver解読の第5回です。）

　絶賛トトリのアトリエプレイ中です！　１周目のクリアが見えましたが、どうやらトゥルーエンドには行かない模様。２周目を見据えて、装備品の錬金に悩む今日この頃です。

それでは本題。今回のキーワードは二重確率行列です。二重確率行列とは何か、どうやって求めるのかについて触れようと思います。

二重確率行列

　前回の話ではBilateral Gridのアフィニティ行列を求めました。論文_*1ではアフィニティ行列を最適化問題に組み込む際、二重確率化という作業を実施しています。二重確率行列を使うと計算式の簡略化を図れるのですが・・・それは次回のお話。今回は二重確率行列の求め方に焦点を当てます。

　二重確率行列（Bistochastic matrix または doubly stochastic matrix）とは何か。Wikiには以下のように書かれています。

各行の和および各列の和がそれぞれ 1 となるような非負の実正方行列
　　　　〜途中略〜
シンクホーンの定理（英語版）では、厳密に正な成分を持つ任意の行列は、適切な対角行列の前方および後方からの乗算によって二重確率行列へと変換することが出来ることが述べられている。

シンクホーンの定理を学べば二重確率行列が計算できるのでしょうか・・・？　結論として、シンクホーンの定理なんぞ知らなくても二重確率行列は計算できます。でもなんかもやもやするので調べてみました。

（参考）*1：The Fast Bilateral Solver

Sinkhorn-Knoppの定理

　調べているとシンクホーン(Sinkhorn)の定理の発展形として、Sinkhorn-Knoppの定理_*2を見つけました。思い切って要約すると、次を主張しています。

台をもつ非負行列は、適切な対角行列の前方および後方からの乗算によって二重確率行列へ変換できる

「厳密に正な行列」から「台をもつ行列」に表現が変わりましたが、台とは聞きなれない単語です。台とは何か、論文を和訳してみました。

（対角成分についての定義）
　 ${\bf A}$ を ${N \times N}$ の大きさもつ行列、 ${\sigma}$ を列 ${(1,2,\cdots ,N)}$ の交換演算子とした時、行列 ${\bf A}$ の要素 ${(a_{1,\sigma(1)}, a_{1,\sigma(2)}, \cdots ,a_{1,\sigma(N)})}$ を ${\sigma}$ に対応した ${\bf A}$ の対角成分と呼ぶ。特に ${\sigma}$ が恒等写像である場合は、主の対角成分と呼ぶ。
（台についての定義）
　以下の条件を満たす場合、 ${\bf A}$ は完全な台（total support）を持つ。
　　・ ${\bf A}$ は非負行列である。（いずれの要素にも負の値をもたない）
　　・ ${\bf A}$ はゼロ行列ではない
　　・全ての非ゼロ要素はいずれかの対角成分に含まれる。
　以下の条件を満たす場合、 ${\bf A}$ は台（support）を持つ。
　　・ ${\bf A}$ は非負行列である。
　　・ ${\bf A}$ は全ての成分が正（0より大きい）となる対角成分をもつ。

　和訳にはレイティング・ランキングの数理―Ｎｏ．１は誰か？を参考にしました。和訳しても意味がわかりません・・・。少なくとも私にはそうです。具体例で考えてみましょう！

$\begin{pmatrix} 0&1&2 \\ 3&4&5 \\ 6&7&8\end{pmatrix}$

主の対角成分は $(a_{1,1}, a_{2,2}, a_{3,3}) = (0, 4, 8)$ のことです。これは当たり前ですよね。では $\sigma$ に対応した対角成分の例を挙げます。 $\sigma=(2,1,3)$ の場合は $(a_{1,2}, a_{2,1}, a_{3,3}) = (1, 3, 8)$ が対角成分になります。 $\sigma=(3,1,2)$ の場合は $(a_{1,3}, a_{2,1}, a_{3,2}) = (2, 3, 7)$ です。他にも３つの対角成分が存在しますが、割愛します。この行列は完全な台を持ちます。

$\begin{pmatrix} 1&2&3 \\ 0&4&5 \\ 0&6&7\end{pmatrix}$

台をもつけど完全ではない例になります。 $\sigma=(1,2,3)$ では $(a_{1,1}, a_{2,2}, a_{3,3}) = (1, 4, 7)$ なので、0を含みません。台をもつ条件を満たします。しかし $\sigma$ をどう選んでも $a_{1,2}=2$ または $a_{1,3}=3$ を含む対角成分は0も含みます。ゆえに完全ではありません。

$\begin{pmatrix} 1&0&0 \\ 0&2&3 \\ 4&0&0\end{pmatrix}$

台をもたない例になります。全ての対角成分が0を含みます。

　以上で台の意味がわかりました。それでは肝心の二重確率行列の求め方に入ります。Sinkhorn-Knopの論文は無視して・・・以下に示すKnightらのアルゴリズムを使います！（要は「台」の意味を知りたかっただけです。）

（参考）*2：Concerning Nonnegative Matrices And Doubly Stochastic Matrices

Knightらによる二重確率化手法

　Knightらが2014年に提案した二重確率行列の計算アルゴリズムを紹介します_*3。行列 $\bf A$ への繰り返し計算により、対角行列 $\bf D,\bf E$ を求めます。 $\bf A'=\bf D \dot \bf A \dot \bf E$ は、 $\bf A$ の対称性を保った二重確率行列となります。

(定理)
$N \times N$ の大きさをもつ行列 $\bf A$ を二重確率化するための行列 $\bf D, \bf E$ を繰り返し計算により求める。
　 $\bf A^{(k)}= (a_{i,j}^{(k)})=\bf D^{(k)} \dot \bf A^{(0)} \dot \bf E^{(k)}$
　 $\bf D^{(k)}= diag(d_1^{(k)}, d_2^{(k)}, \cdots , d_n^{(k)})$ 、ただし $\bf D^{(0)}=\bf I$
　 $\bf E^{(k)}= diag(e_1^{(k)}, e_2^{(k)}, \cdots , e_n^{(k)})$ 、ただし $\bf E^{(0)}=\bf I$
ここで $k$ は繰り返し回数を表わす。また各要素 $d_i,e_j$ は $r_i、c_j$ を用いて以下のように表される。
　 $r_i^{(k)}=||r_i^{(k)}||_1 = \sum^n_{j=1}a_{i,j}^{(k)}$
　 $c_j^{(k)}=||c_j^{(k)}||_1 = \sum^n_{i=1}a_{i,j}^{(k)}$
　 $d_i^{(k)} = 1 / \sqrt{\mathstrut r_i^{(k)}}$
　 $e_j^{(k)} = 1 / \sqrt{\mathstrut c_j^{(k)}}$
この時、以下のことが成り立つ。
　・ $\bf A$ が台をもつ場合、二重確率行列 $\bf S = lim_{k \to \infty}A^{(k)}$ が存在する。
　・ $\bf A$ が完全な台をもつ場合、 $\bf D=lim_{k \to \infty} \bf D^{(k)},\bf E=lim_(k \to \infty) \bf E^{(k)}$ は収束し、二重確率行列 $\bf S=\bf D \dot \bf A \dot \bf E$ が存在する。
　　（補助定理）もしAが対象行列ならば、D=Eとなり、二重確率行列S=DADである。

特に対称行列を扱う場合、繰り返し計算は次のようにシンプルにできます。

$\bf x_{k+1} = \sqrt{\bf x_k / (\bf A \dot \bf x_k})$ 、ただし $\bf x_0=\bf e$ （単位行列）

対称行列は完全な台をもつので、 $\bf x$ は必ず収束します。 $\bf X = diag(\bf x)$ とすると、 $\bf X \dot \bf A \dot \bf X$ は二重確率行列となります。

（参考）*3：A Symmetry Preserving Algorithm For Matrix Scaling