ジョイジョイジョイ

ジョイジョイジョイジョイジョイ

Convolutional LSTM

大学の実験で必要になって実装したのでメモしておきます。

Convolutional LSTM の説明

名前で完全にネタバレしてる感が否めないですが、Convolutional LSTM とは、LSTM の結合を全結合から畳み込みに変更したものです。

例えば画像を RNN に食わすときに、位置情報が失われないので便利です。

動画の次フレームの予測や天気予報などに使えます。

具体的には、以下のような構造になっています。

f:id:joisino:20171027192558j:plain

x は要素ごとの掛け算、 * は畳み込みを表します。

通常の LSTM との差分を赤で書きました。といっても、一部の掛け算が畳み込みになっているだけですが。

peephole の部分だけ要素ごとの掛け算になっていることに注意してください。

実装

Convolutinoal LSTM と、それをビルディングブロックとして使って画像予測のネットワークを実装しました。

github.com

レポジトリです。

python3.5.2 + chainer 3.0.0 で開発しました。

画像予測のネットワークは以下のような構造をしています。([1] より抜粋)

f:id:joisino:20171027151517p:plain

左側のネットワーク(Convolutional LSTM を 2 層並べたもの)に、入力となる画像(例えば 1 フレーム目, 2 フレーム目, 3 フレーム目)を流し込んでいって、データの表現を得ます。

次に、左側の LSTM たちの内部状態を右側の LSTM たちにそれぞれコピーします。

最後に、右側のネットワークにゼロを入力して、出てきた値を concat して、それらを 1x1 で畳み込んで出力を得ることを繰り返します(例えばこれが 4 フレーム目の予測画像, 5 フレーム目の予測画像, 6 フレーム目の予測画像になります)

右側のネットワークの二回目以降の入力は、前回の出力を入れるという実装もありそうです。論文を読んだ限り特に書いてなさそうだったので、今回はゼロを入力するようにしました。

実験

Moving MNIST [2] で実験しました。

計算資源の都合上論文の実装より小さめのネットワークです。

具体的には 2 -> 5x5 Conv -> 128 -> 5x5 Conv -> 64 -> 5x5 Conv -> 64 -> 1x1 Conv -> 2 (右側の forecasting network) という感じです。

以下は生成された画像たちです。

f:id:joisino:20171027145711p:plain

f:id:joisino:20171027143117p:plain

f:id:joisino:20171027143126p:plain

f:id:joisino:20171027143133p:plain

f:id:joisino:20171027143145p:plain

まあいいんじゃないでしょうか。

f:id:joisino:20171027145210p:plain

Loss です。

もう少しちゃんとチューニングする必要がありそうですかね。

だいたい 20 epoch = 8750 iteration = 140000 images の訓練を p2.xlarge (GPU は K80) 上で行って 22 時間くらいでした。

最後に

何か間違いや質問などがあればお願いします。

最後まで読んでいただきありがとうございました。

参考文献

[1] Xingjian Shi, et. al. Convolutional LSTM network: A machine learning approach for precipitation nowcasting. In NIPS, 2015

[2] Moving MNIST

[3] わかるLSTM ~ 最近の動向と共に - Qiita

神絵師がtwitterに上げた神絵を収集する

色んな人がやってそうだけど自分用に作ってみました。

リポジトリです。

github.com

使い方

  • リポジトリをクローンして初期設定します
  • 神絵師のリストを twitter で作って登録します
  • 動かします
  • slack に神絵が流れてきます

動作例

f:id:joisino:20171016205716p:plain

詳しく

動作の流れは以下の通りです。

  1. twitter API を使って、リストから画像を含むツイートを抽出して保存する
  2. 保存した画像から二次元イラストを抽出する
  3. 抽出した画像を slack incoming-webhook をつかってチャンネルに流す

2 が難しいところです。

神絵師が画像を投稿したからといってそれが必ずしも神絵であるとは限りません。

例えば今回の場合、神絵師が食べた焼肉の写真は slack に流れてきてほしくありません。

▲流れてきてほしくないツイートのイメージ(このアカウントは神絵師ではありません)

写真とイラストを区別するのは意外と難しかったです。(何か良いアルゴリズムがあったら教えてください)

ディープラーニングでなんとかしたりできるんでしょうか。

今回は OpenCV の Haar Cascade で nagadomi さんのアニメ顔検出器を使わせてもらいました。

github.com

これでアニメ顔が 1 箇所以上で検出されると、イラストという判定にしました。

体感 accuracy は 7 割くらいですが、 false negative はほとんど無く、slack に流れてくるのは神絵ばかりです。

今回の場合は、網羅的に収集したいというよりは、よさそうな画像をなんとなく眺めたい用途なので、これでよしとします。

今後の課題

ソシャゲガチャ問題

神絵師がソシャゲガチャの結果画像を上げることはよくあります。

今回の場合、ソシャゲガチャの結果画像は slack に流れてきてほしくないですが、アニメ顔が含まれるため防げません。

10 連ガチャの一覧結果はまだ例外でなんとかなりそうですが、特定の一枚の結果を上げている場合、その絵を描いたのも(別の)神絵師なので、弾くのはなかなかむずかしそうです。

最後に

何か間違いや質問などがあればお願いします。

最後まで読んでいただきありがとうございました。

ArtClass(IOI2013)をディープラーニングで解く

はじめに

IOI 2013 オーストラリア大会に Art Class という問題があります。

この問題は、画像データが与えられるのでその画像が

様式1(新造形主義の現代芸術) f:id:joisino:20171003140837p:plain

様式2(印象派の風景画) f:id:joisino:20171003140845p:plain

様式3(表現派のアクション・ペインティング) f:id:joisino:20171003140855p:plain

様式4(カラーフィールド・ペインティング) f:id:joisino:20171003140907p:plain

のいずれであるかを判定する問題です。

正答率が 0.9 以上になると満点が得られます。

IOI にしては珍しい機械学習的な問題であることと、ジャッジが壊れて結果が返ってこなくなったことなどで有名なので、知っている人も多いかもしれません。

問題文やデータは、 http://www.ioinformatics.org/locations/ioi13/contest/ から手に入ります。

普通の解法

例えば 3x365x65 の大きさの窓を作って分散を計算して、それらを使って手で決定木を作るなどすると解けます。

C++ での実装例

gist.github.com

正答率は 99/102 = 0.97... で満点が獲得できました。

ディープラーニングで解く

手で特徴量を作るのは面倒なのでディープラーニングでなんとかして欲しいですよね。

今回は VGG の pretrained model (with ImageNet) を使ってこの問題を解いてみました。(念の為に言っておくと、 IOI の本番では chainer も pretrained model も GPU も使うことはできません)

最終層だけ 4096 x 4 の FC 層に変えて学習しなおしました。

学習データが 4 * 9 = 36 枚しかないのでかなり厳しかったですが、augmentation(crop) してなんとか正答率 92/102 = 0.90... で満点を獲得できました。

適当に作った決定木より精度低い。悲しい。

ImageNet で前学習したモデルを使ったのが筋が悪かったかもしれません。絵画的なデータで前学習したモデルがあれば、もっと楽に満点を達成できたと思います。

テストデータに対しても少しチューニングしたので、汎化性能はさらにもう少し低いと思います。体感 accuracy は 0.8 ~ 0.85 くらいです。IOI 本番でもテストデータに対する正答率は見える(はずだった)のでよしとしましょう。

python + chainer での実装例

環境は python 3.5.2 + chainer 2.1.0 です。

gist.github.com

最後に

正直今回の場合、人力特徴量抽出とハイパーパラメータチューニングの面倒さは同じくらいでした。(本末転倒)

何か間違いや質問などがあればお願いします。

最後まで読んでいただきありがとうございました。

参考文献

画像や問題文は全て http://www.ioinformatics.org/locations/ioi13/contest/ からの引用です。

符号付き整数を二冪で割る

符号付き整数を 2 冪で割るときに右算術シフトを使ってしまって痛い目にあったことのある皆さんこんにちは。僕はあなたの仲間です。

コンピュータで割り算を行うのは一般に重いので、できるだけシフト等の軽い操作で済ませたいですよね。

最近のコンパイラは賢いので適当に書けば最適化してくれますが、何かしらの事情で自分で最適化しなければならいないときのためにやり方を紹介しておきます。

C99 で 32 ビット符号付き整数 x1<<n で割ることを考えます。

符号なし整数なら二冪の割り算 x/(1<<n) と右シフト x>>n は一致するのですが、符号あり整数については、割り算は 0 方向への切り捨てなのに対して、右算術シフトは切り下げとなります。

例えば (-3) / 2 の結果が -1 なのに大して (-3) >> 1 の結果は -2 です。

これは困りました。

簡単な対策として、x が負のときシフトを行う前に (1<<n)-1 を足すという方法があります。

こうすると、負の数のとき切り上げになり、割り算と同じになります。

ですが、if は一般に重いので、できるだけ使わずに済ませたいです。

そこで、全体として以下の手順を踏みます

  1. x31 だけ右算術シフトしたものを y とする。
  2. y32-n だけ右論理シフトしたものを z とする。
  3. xz を足す。
  4. xn 右論理シフトする。

1, 2 が重要です。これによって zx が正のとき 0 に、負のとき (1<<n)-1 になっています。

たとえば n=6 のとき

x = 0xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
y = 00000000000000000000000000000000
z = 00000000000000000000000000000000

x = 1xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
y = 11111111111111111111111111111111
z = 00000000000000000000000001111111

となっています。

正の数のときは最上位ビットが 0 なので、算術シフトによって 0 が引き出され全て 0 で埋まります。論理シフトをしても結果は変わりません。

対して負の数のときは最上位ビットが 1 なので、算術シフトによって 1 が引き出され全て 1 で埋まり、続いて論理シフトによって上位ビット 32-n ビットが全て 0 で埋まります。 1 は全部で n 個並んでいるのでこの値は (1<<n)-1 となります。

例えば以下のように実装できます。

gist.github.com

C99 では論理シフトになるか算術シフトになるかは処理系依存であることに注意してください。

上記のコードは gcc 5.4.0 で動作を確認しています。

おまけ

ちなみに最近のコンパイラでも、ここで紹介した感じになります。

以下は clang 3.8.0x / 64MIPSコンパイルした結果(の一部)です。

sra $2, $1, 31
srl $2, $2, 26
addu $1, $1, $2
sra $1, $1, 6

最後に

何か間違いや質問などがあればお願いします。

最後まで読んでいただきありがとうございました。

参考文献

[1] Henry S. Warren. (2012). Hacker’s Delight (2nd Edition)

シュトラッセンのアルゴリズムとその導出の仕方

 O(n^{3}) よりも小さい時間計算量  O(n^{\log_{2} 7}) \approx O(N^{2.807}) で行列積を計算するシュトラッセンアルゴリズムとその導出方法を紹介します。

シュトラッセンアルゴリズム

 n \times n の行列積  A B = C を考えます。

まず  A, B, C

$$ \begin{pmatrix} A_{11} & A_{12} \\ A_{21} & A_{22} \end{pmatrix} \begin{pmatrix} B_{11} & B_{12} \\ B_{21} & B_{22} \end{pmatrix} = \begin{pmatrix} C_{11} & C_{12} \\ C_{21} & C_{22} \end{pmatrix} $$

というように  \frac{n}{2} \times \frac{n}{2} 行列に分解します。

 n が奇数のときは余った行、列を適当にどちらかに分けてあげてください。

まず  7 つの行列積

$$ \begin{align} P_{1} &= ( A_{11} + A_{22} ) ( B_{11} + B_{22} ) \\ P_{2} &= ( A_{21} + A_{22} ) B_{11} \\ P_{3} &= A_{11} ( B_{12} - B_{22} ) \\ P_{4} &= A_{22} ( B_{21} - B_{11} ) \\ P_{5} &= ( A_{11} + A_{12} ) B_{22} \\ P_{6} &= ( A_{21} - A_{11} ) ( B_{11} + B_{12} ) \\ P_{7} &= ( A_{12} - A_{22} ) ( B_{21} + B_{22} ) \end{align} $$

を計算します。

すると、 C の各ブロックは

$$ \begin{align} C_{11} &= P_{1} + P_{4} - P_{5} + P_{7} \\ C_{12} &= P_{3} + P_{5} \\ C_{21} &= P_{2} + P_{4} \\ C_{22} &= P_{1} + P_{3} - P_{2} + P_{6} \end{align} $$

となります。

この証明は両辺を頑張って展開して一致することを確かめればできます。

ここで、各  P_{i} を計算するときに、再帰的にこのアルゴリズムを適用することを考えます。

各ステップの計算は再帰部分以外は行列の加減算なので、 O(n^{2}) でできます。

再帰は大きさが半分の問題を  7 回呼び出しています。

よって、 Master theorem より、全体の計算量は  O(n^{\log_{2} 7}) \approx O(N^{2.807}) です。

導出方法

 P_{i} C_{ij} の式が与えられれば証明は頑張って計算するだけですが、 P_{i} C_{ij} の式がどうやって出てきたのか分かりません。

ここでは簡単な(とはいっても骨の折れる)導出方法を一つ紹介します。

流れは以下のとおりです。

  1.  2 \times 2 の実行列の掛け算で、要素同士の掛け算を  7 回することを目指す
  2. 平面上の  120 度回転を表す行列  D をとる
  3. 非零ベクトル  u をとる
  4.  u と直交して  u^{\bot} D u = 1 を満たす横ベクトル  u^{\bot} をとる
  5.  X = u u^{\bot} とする
  6.  \{ D, X, D^{-1}XD, DXD^{-1} \} \{ D^{-1}, X, D^{-1}XD, DXD^{-1} \} 2 \times 2 の行列の基底になっている
  7.  A, B をそれぞれ上の基底で表現してかけると  7 回ですむ
  8.  D は実は固有多項式が同じなら同じ議論が成り立つので、簡単な行列をとってみる
  9.  u も簡単なベクトルをとってみる
  10. すると実は 7 はシュトラッセンアルゴリズムと同じだった

長いですね。

一つずつ見ていきましょう。

回転行列  D

 120 度 ( 2/3 \pi) の回転行列

$$ D = \begin{pmatrix} -\frac{1}{2} & -\frac{\sqrt{3}}{2} \\ \frac{\sqrt{3}}{2} & -\frac{1}{2} \end{pmatrix} $$

を取ります。

 tr(D) = -1,  det(D) = 1 なので、固有多項式 \lambda^{2} + \lambda + 1 です。

ケーリーハミルトンの定理より  D^{2} + D + E = 0 \cdots (1)が成り立ちます。

両辺に  D をかけると  D^{3} + D^{2} + D = 0 \cdots (2) が成り立ちます。

(1) と (2) の左辺どうしは等しいので  D^{3} = E \cdots (3) となります。これは、 D 120 度の回転行列であることを考えると当たり前のことです。

(3) より  D^{-1} = D^{2} \cdots (4) です。

よって、 (1) は  D^{-1} + D + E = 0 \cdots (5) と書き換えられます。

(5) より  tr(D^{-1}) = -tr(D) - tr(E) = 1 - 2 = -1 \cdots (6) となります。

非零ベクトル  u

非零ベクトル  u をとります。

 u と直交する非零の横ベクトル  u' をとります。

 D の固有多項式の実根は存在しないので、 u D固有ベクトルではありません。

よって、 Du u と別方向を向いているはずなので、  u'Du \neq 0 です。

 u^{\bot} = \frac{u'}{u'Du} とします。

 u^{\bot}u と直交し  u^{\bot} D u = 1 となります。

行列  X

 X = u u^{\bot} をとります。ここで、

 X^{2} = (u u^{\bot}) (u u^{\bot}) = u (u^{\bot} u) u^{\bot} = 0 \cdots(7)

 XDX = (u u^{\bot}) D (u u^{\bot}) = u (u^{\bot} D u) u^{\bot} = u u^{\bot} = X \cdots(8)

 XD^{-1}X = (u u^{\bot}) D^{-1} (u u^{\bot}) = u (u^{\bot} D^{-1} u) u^{\bot} = -u u^{\bot} = -X \cdots(9)

が成り立ちます。

一次独立性

前々節でも言ったように、 D の固有多項式の実根は存在しないので、 u D固有ベクトルではありません。

よって、 u Du は一次独立です。

同様に  u^{\bot} u^{\bot} D は一次独立です。

よって、これらの積  u \cdot u^{\bot}, u \cdot u^{\bot} D, Du \cdot u^{\bot}, Du \cdot u^{\bot} D 2 \times 2 行列の基底になります。

これらに両側から  D を掛けます。

 D^{-1} = D^{2} \cdots (4) を考えると、  DXD, DXD^{-1}, D^{-1}XD, D^{-1}XD^{-1} となります。

 D は正則なので、これらも  2 \times 2 行列の基底になります。

これらを足し合わせると  (D + D^{-1}) X (D + D^{-1}) = (-E) X (-E) = X となります。ひとつ目の等号は  D^{-1} + D + E = 0 \cdots (5) から従います。

よって、基底のうち一つを  X で置き換えても基底になります。

これより、  X, DXD^{-1}, D^{-1}XD は一次独立です。

基底表現

 X, DXD^{-1}, D^{-1}XD のトレースはそれぞれ、

 tr(X) = tr(u^{\bot}u) = tr(uu^{\bot}) = 0

 tr(DXD^{-1}) = tr(XD^{-1}D) = tr(X) = 0

 tr(D^{-1}XD) = tr(XDD^{-1}) = tr(X) = 0

となります。

よって、 \{ X, DXD^{-1}, D^{-1}XD \} から生成される  2 \times 2 行列の部分空間の要素のトレースは  0 です。

 tr(D) = -1, tr(D^{-1}) = -1 より、これらはそれぞれ  X, DXD^{-1}, D^{-1}XD と一次独立です。

よって、 \{ D, X, D^{-1}XD, DXD^{-1} \} \{ D^{-1}, X, D^{-1}XD, DXD^{-1} \} はそれぞれ  2 \times 2 行列の基底になっています。

$$ A = a_{1} D + a_{2} X + a_{3} D^{-1}XD + a_{4} DXD^{-1} $$

$$ B = b_{1} D^{-1} + b_{2} X + b_{3} D^{-1}XD + b_{4} DXD^{-1} $$

と表します。このとき基底同士の掛け算は

f:id:joisino:20170903160923p:plain

となります。ここでは (7), (8), (9) の結果をたくさん使っています。(表は [2] より引用)

同じ結果が同列、同行に現れているので、これらはまとめられます。

つまり、

$$ \begin{align} AB &=& a_{1} &\times& b_{1} &\times& id \\ &+& a_{2} &\times& ( b_{1} + b_{4}) &\times& XD^{-1} \\ &+& a_{3} &\times& ( b_{1} + b_{2}) &\times& D^{-1}X \\ &+& a_{4} &\times& ( b_{1} + b_{3}) &\times& DXD \\ &+& ( a_{1} - a_{4}) &\times& b_{2} &\times& DX \\ &+& ( a_{1} - a_{2}) &\times& b_{3} &\times& XD \\ &+& ( a_{1} - a_{3}) &\times& b_{4} &\times& D^{-1}XD^{-1} \cdots (10) \end{align} $$

となります。

もう少し簡単に

ここまで  D の性質としては  tr(D) = -1,  det(D) = 1 であることしか使っていません。よって、たとえば

$$ D = \begin{pmatrix} 0 & -1 \\ 1 & -1 \end{pmatrix} $$

ととっても構いません。このとき

$$ D^{-1} = \begin{pmatrix} -1 & 1 \\ -1 & 0 \end{pmatrix} $$

です。

また、 u にも任意性がありましたが、例えば  u = ( 1, 0 )^{{\rm T}} ととります。

このとき

$$ X = \begin{pmatrix} 0 & 1 \\ 0 & 0 \end{pmatrix} $$

となります。

このようにとることで上記の議論を  2 \times 2 行列の要素がブロック行列でも適用できることになります。( 1単位行列で、 0 は零行列です)

また、計算すると

$$ D^{-1}XD = \begin{pmatrix} -1 & 1 \\ -1 & 1 \end{pmatrix} $$

$$ DXD^{-1} = \begin{pmatrix} 0 & 0 \\ -1 & 0 \end{pmatrix} $$

となります。 A, B の要素ごとに考えると、

 X, DXD^{-1}, D^{-1}XD はトレースが  0 なので、  a_{1} = - tr(A) = - a_{11} - a_{22} となります。

また、 X, D^{-1}XD (2,2) 要素が  0 なので、 a_{3} = a_{22} + a_{1} = - a_{11} となります。

残りは、それぞれ異なる場所に  1, -1 が立っているだけなので、それぞれ

 a_{2} = a_{12} + a_{1} - a_{3} = a_{12} - a_{22}

 a_{4} = -( a_{21} - a_{1} + a_{3}) = - a_{21} - a_{22}

となります。

同様に

 b_{1} = -b_{11} - b_{22}

 b_{2} = b_{11} + b_{12}

 b_{3} = b_{22}

 b_{4} = b_{11} - b_{21}

となります。

また、瑣末な計算なのでここでは省きますが、 id, XD^{-1}, D^{-1}X, DXD, DX, XD, D^{-1}XD^{-1} もそれぞれ  1, 0, -1 からなる簡単な  2 \times 2 行列になっています。

最後に (10) で表したようにそれぞれ係数を掛けて足し合わせると、シュトラッセンアルゴリズムが出てくるというわけです。

最後に

何か間違いや質問などがあればお願いします。

もっと簡単な導出があれば教えていただけると嬉しいです。

最後まで読んでいただきありがとうございました。

参考文献

[1] Volker Strassen. (1969) Gaussian elimination is not optimal.

[2] Christian Ikenmeyer. Vladimir Lysikov. (2017) Strassen’s 2x2 matrix multiplication algorithm: A conceptual perspective.

[3] Master theorem - Wikipedia

ケイリーの公式の証明6種類

ケイリーの公式の証明たちの紹介です。

ケイリーの公式とは

ケイリーの公式とは  n 頂点のラベル付きの木の総数  T(n) T(n) = n^{n-2} であるという公式のことです。

ここで、ラベル付きであるとは、それぞれの頂点を区別するということです。

たとえば  n = 3 のとき、頂点を区別しない場合は長さ  3 のパスのみの  1 通りですが、ラベル付きの木の場合は  1 - 2 - 3,  2 - 1 - 3,  1 - 3 - 2 3 (= 3^{3-2} = T(3)) 通りです。

証明 1 (プリューファーコード) [1]

おそらく一番有名な証明です。

 n 頂点のラベル付きの木の集合から  \{ 1, 2, \ldots, n \}^{n-2} への全単射を以下のように構成します。

最もラベルが小さい葉を木から取り除き、その葉と繋がっていた頂点のラベル  a_{1} を数列の最初の値とします。

続けて、最もラベルが小さい葉を木から取り除き、その葉と繋がっていた頂点のラベル  a_{2} を数列の  2 番目の値とします。

以下同様に頂点が  2 つになるまで操作を続けます。

こうしてできた数列  a_{1}, a_{2}, \ldots, a_{n-2} が木の値となります。

この数列をプリューファーコードといいます。

例えば、以下の木のプリューファーコードは  3, 2, 6, 6 となります。

f:id:joisino:20170819151034p:plain

これの逆写像は以下のように与えられます。

まず、辺のない  n 頂点を用意します。

次に、  a_{1}, a_{2}, \ldots, a_{n-2} に含まれない値のうち、最も小さいものを  b_{1} とし、  a_{1} b_{1} を繋ぎます。

続けて、  b_{1}, a_{2}, a_{3} \ldots, a_{n-2} に含まれない値のうち、最も小さいものを  b_{2} とし、  a_{2} b_{2} を繋ぎます。

続けて、  b_{1}, b_{2}, a_{3}, a_{4} \ldots, a_{n-2} に含まれない値のうち、最も小さいものを  b_{3} とし、  a_{3} b_{3} を繋ぎます。

以下同様に木を構成していきます。

最後に、 b_{1}, b_{2}, \ldots, b_{n-2} に含まれていない  2 頂点を繋いで完成です。

これが逆写像になっていることは以下のようにわかります。

まず、プリューファーコードに登場する頂点は元の木の内点で、かつ元の木の内点を全て尽くしています。

なぜなら

  •  a_{i} i 番目の操作の時点で葉とつながっている点なので、これが葉であると、操作後に頂点数が 1 の成分になってしまうので矛盾です。よって全ての  a_{i} は内点です。

  • 全ての内点は、操作のどこかの時点で葉になり、 2 頂点以外は取り除かれていきます。少なくともその頂点が葉になった瞬間においてコード中に登場するはずなので、全ての内点を尽くしています。

以上より、逆に  a_{1}, a_{2}, \ldots, a_{n-2} に含まれない値は初期状態において葉であって、全ての葉を尽くしています。

そのうち最もラベルの小さい頂点は順写像において最初に取り除かれた頂点で、 a_{1} と繋がっていたはずです。

以下帰納的に順写像の操作に対応する辺が見つかっていきます。

最後は、順写像の操作で取り除かれなかった頂点を結べば完成となります。

全単射が存在するということは要素の数が等しいということです。

 \{ 1, 2, \ldots, n \}^{n-2} の要素数はもちろん  n^{n-2} なので  T(n) = n^{n-2} となります。

証明 2 (もうひとつの全単射) [2]

もうひとつ全単射を使った証明を紹介します。

今度は  \{ 1, 2, \ldots, n \}^{n} から、  n 頂点からなるラベル付きの木であって  1 頂点に赤い印がついていて  1 頂点が星形である木への全単射を構成します。

ここで、赤の頂点と星の頂点は同じ頂点でもかまいません。

ラベル付きの木を固定すると赤の頂点と星の頂点の選び方は  n^{2} 通りあるので、値域の要素数 T(n) n^{2} で、定義域の要素数 n^{n} なので、これらが等しければケイリーの公式の証明が完了します。

数列  (a_{1}, a_{2}, \ldots, a_{n}) を入力とします。

 1 \le i \le n について  i から  a_{i} へ辺があるような  n 頂点  n 辺からなるグラフを考えます。

全ての出次数が  1 であるグラフは functional graph と呼ばれ、サイクルに木がついた連結成分がいくつかあるような形になります。

例えば  7, 4, 8, 1, 7, 4, 4, 3 のときグラフは以下のようになります。

f:id:joisino:20170819161216p:plain

ここで、サイクル上の頂点をソートして並べたものを  b_{1}, b_{2}, \ldots, b_{k} とします。

上の例の場合  1, 3, 4, 7, 8 となります。

出力の木としては、  a_{b_{1}}, a_{b_{2}}, \ldots, a_{b_{k}} をこの順に並べてパスとし、残りは functional graph と同様に繋げます。赤の頂点は  a_{b_{1}} を選び、星の頂点は  a_{b_{k}} を選びます。

上の例の場合は以下のようになります。

f:id:joisino:20170819162001p:plain

写像は以下のように構成します。

まず、木のうち、赤の頂点から星の頂点までの(一通りしかない)パスをとって、 c_{1}, c_{2}, \ldots, c_{k} と順に並べます。

パス上の頂点のラベルをソートした数列を  b_{1}, b_{2}, \ldots, b_{k}とします。

 1 \le i \le k について、元の数列の  b_{i} 番目の値が  c_{i} であることが分かります。

パス以外の頂点は、パス上の頂点を根とした根付き木を考え、頂点  v の親が  u のとき、元の数列の  v 番目の値が  u となることが分かります。

これで復元できました。

証明 3 (double counting) [2] [3]

この証明もかなり有名な証明です。

wikipedia にもこの証明が載っています。

頂点にも辺にもラベルが付いた  n 頂点からなる根付き木の総数  U(n) を二通りの方法で数えます。

一つは、 T(n) を使う方法です。ラベル付きの木を固定したとき、根の選び方が  n 通り、辺へのラベルの付け方が  (n-1)! 通りあるので、 U(n) = T(n) \cdot n \cdot (n-1)! となります。

もうひとつは辺を  1 つずつ追加していって逐次的に木を構成していく方法です。 i 番目に追加した辺のラベルを  i とします。

辺のない  n 頂点からなるグラフを用意します。

最初の辺は、親の選び方が  n 通りあって、子の選び方は自分以外の  n-1 通りです。

 2 つ目の辺は、親の選び方がこれも  n 通りあります。子の選び方は、根でない頂点を選ぶとその頂点の親が  2 つになってしまうのでだめです。自分が属している木の根も、ループができてしまうのでだめです。よって、子の候補は  n-2 個となります。

同様に、 1 \le i \le n-1 について  i 番目の辺は、親の選び方が  n 通り、子の選び方が  n-i 通りとなります。

よって、 U(n) = \Pi_{i=1}^{n-1} n (n-i) = n^{n-1} \cdot (n-1)! となります。

これが  T(n) \cdot n \cdot (n-1)! と等しいので  T(n) = n^{n-2} となります。

証明 4 (行列木定理) [2]

行列木定理 を使った証明です。

ラベル付きの木の総数と  n 頂点完全グラフ  K_{n} の全域木の総数は同じです。

 K_{n}ラプラシアン行列は

$$ \displaystyle \left[ \begin{array}{rrr} n-1 & -1 & -1 & \ldots & -1 \\ -1 & n-1 & -1 & \ldots & -1 \\ -1 & -1 & n-1 & \ldots & -1 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ -1 & -1 & -1 & \ldots & n-1 \end{array} \right] $$

となります。行列木定理より、 T(n) はこの行列の任意の余因子と一致します。

 (1, 1) 余因子をとりましょう。

 (1, 1) 余因子はもとの行列から  1 行目と  1 列目を取り除いた  (n-1) \times (n-1) 行列の行列式です。

まず、 2 行目から  n-1 行目を  1 行目に足します。すると  1 行目が全て  1 になります。

$$ \displaystyle \left| \begin{array}{rrr} 1 & 1 & 1 & \ldots & 1 \\ -1 & n-1 & -1 & \ldots & -1 \\ -1 & -1 & n-1 & \ldots & -1 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ -1 & -1 & -1 & \ldots & n-1 \end{array} \right| $$

次に、 1 行目を  2 行目から  n-1 行目に足します。

$$ \displaystyle \left| \begin{array}{rrr} 1 & 1 & 1 & \ldots & 1 \\ 0 & n & 0 & \ldots & 0 \\ 0 & 0 & n & \ldots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \ldots & n \end{array} \right| $$

この式の値は  n^{n-2} なので、  T(n) = n^{n-2} となります。

証明 5 (漸化式) [2]

より一般的な問題を漸化式と帰納法によって解きます。

ラベルの集合の、要素数  k の部分集合  A を任意にとります。

 A に含まれている頂点が互いに異なる連結成分に属するようなラベル付き森の数を  S(n,k) とします。

これは  A の要素数にしかよらないので well-defined です。

便宜上  S(0,0) = 1, S(n,0) = 0 ( n > 0 ) とします。

 S(n,k) = k \cdot n^{n-k-1} (*) であることを示します。

 T(n) = S(n,1) なので、これが示せればケイリーの公式は示されます。

ラベルの集合を  \{ 1, 2, \ldots, n \} とし、 A = \{ 1, 2, \ldots, k \} としたときに  S(n, k) で数える森を  1 つ固定します。

この森で頂点  1 と繋がっている頂点を  B = \{ a_{1}, a_{2}, \ldots, a_{l} \} とします。

このときこの森から頂点  1 を取り除いてできる森は、ラベルの集合を  \{ 2, 3, \ldots, n \} とし  A = \{ 2, 3, \ldots, k, a_{1}, a_{2}, \ldots, a_{l} \} としたときに  S(n-1,k+l-1) で数える要素になります。

逆に ラベルの集合を  \{ 2, 3, \ldots, n \} とし  A = \{ 2, 3, \ldots, k, a_{1}, a_{2}, \ldots, a_{l} \} としたとき、新たに頂点  1 を加えて  a_{1}, a_{2}, \ldots, a_{l} とつなげると、先ほど考えていた  S(n, k) で数えた森に対応します。

つまり、 S(n,k) で数える森は、 A と交わらない頂点集合  B ごとに  S(n-1,k+|B|-1) 個の森と対応づいているということです。

以上の考察から

 \displaystyle S(n,k) = \sum_{l=0}^{n-k} \binom{n-k}{l} S(n-1, k+l-1)

となります。これは  B として、 A 以外の要素 n-k 個から  l 個選ぶというのを  l = 0, 1, \cdots n-k について足しあわせているということです。

この漸化式を使って帰納法により (*) を示しましょう。

 
\begin{align*}
S(n,k) &= \sum_{l=0}^{n-k} \binom{n-k}{l} S(n-1, k+l-1) & \\
&= \sum_{l=0}^{n-k} \binom{n-k}{l} (k+l-1) (n-1)^{n-k-l-1} & ({\rm 帰納法の仮定より})\\
&= \sum_{l=0}^{n-k} \binom{n-k}{l} (n-1-l) (n-1)^{l-1} & (l = n - k - l {\rm と変換}) \\
&= \sum_{l=0}^{n-k} \binom{n-k}{l} (n-1) \cdot (n-1)^{l-1} - \sum_{l=0}^{n-k} \binom{n-k}{l} l \cdot (n-1)^{l-1} & ({\rm 分配}) \\
&= \sum_{l=0}^{n-k} \binom{n-k}{l} (n-1)^{l} - \sum_{l=0}^{n-k} \binom{n-k}{l} l \cdot (n-1)^{l-1} & \\
&= (n-1+1)^{n-k} - \sum_{l=0}^{n-k} \binom{n-k}{l} l \cdot (n-1)^{l-1} & ({\rm 二項定理})\\
&= n^{n-k} - \sum_{l=0}^{n-k} \binom{n-k}{l} l \cdot (n-1)^{l-1} & \\
&= n^{n-k} - \sum_{l=1}^{n-k} \binom{n-k}{l} l \cdot (n-1)^{l-1} & (l=0 {\rm のとき} 0) \\
&= n^{n-k} - (n-k) \sum_{l=1}^{n-k} \binom{n-k-1}{l-1} \cdot (n-1)^{l-1} & ({\rm 二項係数をほぐす})\\
&= n^{n-k} - (n-k) \sum_{l=0}^{n-k-1} \binom{n-k-1}{l} \cdot (n-1)^{l} & (l=l+1 {\rm と変換})\\
&= n^{n-k} - (n-k) (n-1+1)^{n-k-1} & ({\rm 二項定理})\\
&= n^{n-k} - (n-k) n^{n-k-1} & \\
&= k \cdot n^{n-k-1} &
\end{align*}

証明 6 (母関数とラグランジュの反転公式) [1]

 n 頂点からなるラベル付きの根付き木の総数を  R(n) とします。便宜上  R(0) = 0 とします。

根なし木を固定したとき根の選び方は  n 通りなので、 R(n) = n T(n) となります。 R(n) = n^{n-1} を示すのが目標です。

根の選び方は  n 通りあります。

根の子の個数を  m に固定します。

子に順番があると考え、それぞれの部分木の大きさが左から順に  k_{1}, k_{2}, \ldots, k_{m} (k_{1} + k_{2} + \ldots + k_{m} = n-1) であるとします。

根の子の部分木について、どのラベルを割り当てるかが  \frac{(n-1)!}{k_{1}! k_{2}! \ldots k_{m}!} 通りあります。

各部分木について根付き木の作り方はそれぞれ  R(k_{1}), R(k_{2}), \ldots, R(k_{m}) 通りあります。

子の順番を取り除くために、これらの積を  m! で割る必要があります。

以上より

 
\begin{align*}
R(n)
&= n \sum_{m=0} \frac{1}{m!} \sum_{k_{1} + k_{2} + \ldots + k_{m} = n-1} \frac{(n-1)!}{k_{1}! k_{2}! \ldots k_{m}!} R(k_{1}) \cdot R(k_{2}) \cdot \ldots \cdot R(k_{m}) \\
&= \sum_{m=0} \frac{1}{m!} \sum_{k_{1} + k_{2} + \ldots + k_{m} = n-1} \frac{n!}{k_{1}! k_{2}! \ldots k_{m}!} R(k_{1}) \cdot R(k_{2}) \cdot \ldots \cdot R(k_{m}) \\
\end{align*}

という漸化式が成り立つことが分かります。

次に  R(n) の指数型母関数

 \displaystyle f(x) = \sum_{n=0}^{\infty} \frac{ R(n) x^n }{ n! }

を考えます。

 
\begin{align*}
x e^{f(x)}
&= x \sum_{m=0}^{\infty} \frac{ f(x)^m }{ m! } \\
&= x \sum_{m=0}^{\infty} \frac{1}{m!} (\sum_{n=0}^{\infty} \frac{ R(n) x^n }{ n! })^m \\
&= x \sum_{m=0}^{\infty} \frac{1}{m!} \sum_{k_{1}=0}^{\infty} \sum_{k_{2}=0}^{\infty} \cdots \sum_{k_{m}=0}^{\infty} \frac{1}{k_{1}! k_{2}! \ldots k_{m}!} R(k_{1}) \cdot R(k_{2}) \cdot \ldots \cdot R(k_{m}) x^{k_{1} + k_{2} + \ldots + k_{m}} \\
&= \sum_{m=0}^{\infty} \frac{1}{m!} \sum_{k_{1}=0}^{\infty} \sum_{k_{2}=0}^{\infty} \cdots \sum_{k_{m}=0}^{\infty} \frac{1}{k_{1}! k_{2}! \ldots k_{m}!} R(k_{1}) \cdot R(k_{2}) \cdot \ldots \cdot R(k_{m}) x^{k_{1} + k_{2} + \ldots + k_{m}+1} \\
&= \sum_{m=0}^{\infty} \frac{1}{m!} \sum_{k_{1}=0}^{\infty} \sum_{k_{2}=0}^{\infty} \cdots \sum_{k_{m}=0}^{\infty} \frac{(k_{1} + k_{2} + \ldots + k_{m} + 1)!}{k_{1}! k_{2}! \ldots k_{m}!} R(k_{1}) \cdot R(k_{2}) \cdot \ldots \cdot R(k_{m}) \frac{x^{k_{1} + k_{2} + \ldots + k_{m}+1}}{(k_{1} + k_{2} + \ldots + k_{m} + 1)!} \\
&= f(x)
\end{align*}

最後の等式は先ほどの漸化式から従います。

よって、 R(n) の指数型母関数  f(x) f(x) = x e^{f(x)} という単純な方程式に従うことがわかりました。

次にラグランジュの反転公式を使います。

ラグランジュの反転公式とは、べき級数

 \displaystyle f(x) = a + \frac{a_{1}}{1!} x + \frac{a_{2}}{2!} x^{2} + \ldots

 \displaystyle f(x) = a + x \phi( f(x) )

という方程式を満たすとき、

 \displaystyle a_{n} = \left| \frac{{\rm d}^{n-1}}{{\rm dt}^{n-1}} (\phi(t)^n) \right|_{t=a}

を満たすというものです。

これに先ほどの方程式を当てはめます。

(本当は収束性などを調べないといけないはずですが、参考文献では特に言及されていませんでした。実際  R(n) = n^{n-1} なので収束半径は  \frac{1}{e} になるはずなのですが、答えを使わずに証明する良い方法が思いつかなかったのでここでは省略します。僕が勘違いしていたり、良い方法があったりしたら教えてください)

ここで、 \phi(t) = e^{t}, a = 0 なので、

 \displaystyle \frac{{\rm d}^{n-1}}{{\rm dt}^{n-1}} e^{nt} = n^{n-1} e^{nt}

よって、 t = a = 0 を代入すると  R(n) = n^{n-1} となり示されました。

最後に

何か間違いや質問などがあればお願いします。

他の証明もあれば教えていただけると嬉しいです。

最後まで読んでいただきありがとうございました。

参考文献

[1] パラモノヴァ(著). ヴィンベルク(著). コーハシ(著). 武部尚志(訳). (2007). モスクワの数学ひろば 3 代数篇/対称性・数え上げ.

[2] Aigner Martin. Ziegler Günter M. (2013). Proofs from THE BOOK 4th edition. pp.201-206

[3] ケイリーの公式 - Wikipedia

前処理O(n)クエリO(1)のLCAと静的RMQ

時間計算量 <O(n), O(1)> の LCA(Lowest Common Ancestor) と RMQ(Range Minimum Query) を C++ で実装しました。

アルゴリズムの解説はDさんのスライド [1] LCA and RMQ ~簡潔もあるよ!~ がとても分かりやすいのでそちらを参照してください。

概要だけ説明します。

LCA の概要

LCA は頂点を dfs 順で訪れた順に並べると深さの列の RMQ に帰着されます。このことは [2] 蟻本 などに載っています。

この列は隣り合う数の差がちょうど  1 になっています。

この列を  \frac{\log n}{2} 個ずつのブロックに分け、それぞれのブロック内の最小値を求めます。

ブロックの数は  \frac{2n}{\log n} 個になるので、ブロックの区間の最小値を求めるクエリは sparse table を使うと前処理  \frac{2n}{\log n} \log \frac{2n}{\log n} = O(n)、クエリ  O(1) で処理できます。

ブロックの中についてですが、各ブロックは 1 つ目の数がわかればあとは +1 されるか -1 されるかの 2 通りです。この遷移の種類は高々  2^{\frac{\log n}{2} - 1} = O(\sqrt{n}) 個しかありません。

この全ての種類について、最初の数が  0 だと思って愚直に前計算すると全て合わせて  O(\log^{2} n \cdot \sqrt{n}) で計算できます。

あとはクエリが来た時に、完全に含まれるブロックを sparse table で処理して、はみ出た部分を前計算した表を参照して答えればよいです。

RMQ の概要

列を Cartesian Tree という木に変換します。

Cartesian Tree は、列のインデックスをノードの値にもつ二分探索木で、列の値についてヒープ条件を満たしています。ヒープ条件を満たしているとは、列の値について、親の値が子の値以下であるということです。

Treap のような感じです。

すると、RMQ はこの Cartesian Tree 上の LCA クエリに帰着されることがわかると思います。

Cartesian Tree の構成の仕方は [1] には載っていないので少し詳しく説明します。

Cartesian Tree を O(n) で構成する

列を左から右に見ていって、順に Cartesian Tree を構成していきます。

今まで構成した Cartesian Tree の頂点の親を配列で管理し、根から右の子を辿るパスをスタックで管理します。

新しく列の最後に値を追加するとします。

このとき、二分探索木の条件から、新しい頂点はスタックで管理されているパス上の頂点の右の子になるか、根になるかのどちらかです。

どの頂点の子になるかは、スタックの頂点の値よりも新しい頂点の値の方が小さければ、スタックから頂点を取り除くことを繰り返し、最終的にスタックの一番上に残っている頂点が新しい頂点の親になります。スタックが空ならば新しい頂点は根になります。

新しい頂点のせいで親から引き離された頂点は、新しい頂点の左の子にすればよいです。

例えば 4, 2, 5, 1, 3, 8, 7, 9 の後ろに 6 を付け加える場合は以下のようになります。

ノードは (index, val) のように表しています。赤いノードがスタックで管理されているパスです。

f:id:joisino:20170813142841p:plain

f:id:joisino:20170813142849p:plain

実装

C++ で実装しました。

LCA

gist.github.com

RMQ

gist.github.com

LCA は GRL_5_C に投げたところ、 10^{5} 頂点  10^{5} クエリで AOJ 上で 0.06s でした。入出力にもある程度時間がかかっていると思います。

RMQ の方は、DSL_3_D (スライド最小値) に投げたところ、  10^{5} 要素  10^{5} クエリで AOJ 上で 0.34s でした。 RMQ は何段階も帰着させている分線形とはいえ遅いです。

メルセンヌツイスタで生成した乱数に対して RMQ を実行すると、手元の環境だと

  •  10^{6} 要素  10^{6} クエリで 0.425s
  •  10^{7} 要素  10^{7} クエリで 5.206s

でした。

最後に

何か間違いや質問などがあればお願いします。

最後まで読んでいただきありがとうございました。

参考文献

[1] LCA and RMQ ~簡潔もあるよ!~

[2] 秋葉拓哉. 岩田陽一. 北川宜稔. (2012). プログラミングコンテストチャレンジブック [第2版]. pp.294 - 295.