ジョイジョイジョイ

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

Dilated Convolution

Dilated Convolution を chainer で実装しました。

Dilated Convolution の説明

Dilated Convolution は、フィルターとの積を取る相手の間隔をあける畳み込みのことです。

例えば、以下のような画像において、 12 を中心に 3 x 3 の普通の畳み込みフィルターを適用すると、 6, 7, 8, 11, 12, 13, 16, 17, 18 との積を取って和を取ると思います。

0 1 2 3 4
5 6 7 8 9
10 11 12 13 14
15 16 17 18 19
20 21 22 23 24

3 x 3 の dilate = 2 の Dilated Convolution フィルターを 12 を中心に適用すると、0, 2, 4, 10, 12, 14, 20, 22, 24 と 1 つおきに取ってきて、それらに 3 x 3 の畳み込みフィルターを適用します。

これは 5 x 5 の通常の畳み込みフィルターのうち、9 箇所以外を 0 に固定したものと見ることもできます。

dilate = 1 の Dilated Convolution は普通の畳込みになります。

dilate が 1 増えるごとに間隔が 1 ずつ開きます。

縦横で異なる dilate を使用する場合もあります。

この Dilated Convolution の何が嬉しいかというと、受容野(あるマスに影響する入力部分)を簡単に増やすことができることです。

例えば、 3 x 3 の畳み込み層を 10 層適用すると、出力の 1 つのマスにある情報は、最初の層の 21 x 21 の部分の情報をまとめたものになります。

入力画像がとても大きい場合、これではごく局所的な情報しか集約できていません。

ここで、簡単に受容野を増やすには

  1. 層の数を増やす
  2. フィルターを大きくする
  3. プーリング層を使う

などの方法があります。

1, 2 については、受容野の一辺の長さは層数とフィルターの一辺の長さに対して線形なので、画像がとても大きい場合はあまり有用ではありません。

3 は簡単に情報を集約できるので有用で、昔からよく使われてきました。

ただ、3 の欠点の一つとして、プーリング層を入れると画像がその分小さくなってしまうことがあります。

例えば画像を入力として、そのうち人間が写っている部分を出力したいといったタスクを考えると、できるだけ元の入力と同じ大きさの画像を出力したいです。

また、ネットワークを深くしたいときに画像の大きさによって限界がきてしまうのも困ります。

そこで、Dilated Convolution を使うと、この欠点なしに情報を集約できます。

よく使われるのは、Dilate を 1, 2, 4, 8, 16, …, 512 と指数的に増やす方法です。

こうすると、パラメータ数は層数に線形なのに対して、受容野の大きさを指数的に増やすことができ、画像もそこまで小さくなりません(パディングで対応できる程度です)。

実際に、 Dilated Convolution は大域的な情報が必要なタスクでは有用なようです。

実装

gist.github.com

chainer で実装してみました。

最初は愚直な 7 重ループで実装したのですが、重すぎて学習が始まらなかったので書き直しました。

numpy 配列をまとめて計算するのとランダムアクセスするのが同じ計算時間でできると思うと痛い目をみるようです。

このへんの細かい仕様はよく分かってないのでまた調べたいです。

二度目の実装にあたっては chainer の dilated_convolution_2d を大いに参考にしました。

相違点として気をつけるべきなのが、変形後の x の次元が、chainer の方では (batch, in_layer, conv_y, conv_x, out_height, out_width) なのに対し、この実装は (batch, in_layer, out_height, out_width, conv_y, conv_x) となってることです。

こちらの方が直感的だとは思ったのですが、よく考えると画像の大きさよりもフィルターの大きさの方が小さいので、 chainer の実装方針でやったほうが効率は良いはずです。

他にも要因はあるかと思いますが、実際に速度は 4 倍ほど chainer の実装の方が速いです。

また、パラメータの初期値でもハマりました。

最初は全てのパラメータを平均 0 分散 1 の正規分布で初期化していたのですが、さすがに大きすぎたようでほとんど精度が出ませんでした。

chainer 標準の initializer を使うとうまくいきました。

initializer も中身をよく理解していないのでまた調べたいです。

実験

MNIST の分類を行いました。

Dilated Convolution

f:id:joisino:20170712164414p:plain

通常の Convolution

f:id:joisino:20170712164439p:plain

ほとんど変わりませんね。

むしろ最終的に dilate = 1 の方が精度が出ていて残念です。

MNIST の分類程度では違いは出ないということなのでしょうか。

Dilated Convolution のほうでもちゃんと学習が進んでくれたのでよしとします。

最後に

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

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

参考文献

[1] F. Yu et al. (2015). Multi-Scale Context Aggregation by Dilated Convolutions

[2] Dilated Convolutions and Kronecker Factored Convolutions ( http://www.inference.vc/dilated-convolutions-and-kronecker-factorisation/ )

[3] van den Oord et al. (2016). WaveNet: A Generative Model for Raw Audio

Batch Normalization

Batch Normalization [1] を chainer で実装しました。

色々な場合に適用できて、学習速度が速くなったり汎化性能が上がったりするすごいテクです。

Batch Normalization の説明

上層(出力層に近い層)の入力は、当然下層(入力層に近い側)のパラメータに依存します。

学習が進むにしたがって下層のパラーメータは変化するので、それにしたがって上層の入力の分布が変化します。

このような中間層の入力の分布の変化を Internal Covariate Shift と呼びます

入力の分布が変化した層はそれに合わせてパラメータも学習しなおさなければなりません。

これが、学習の速度を下げる要因になります。

Batch Normalization は、Internal Covariate Shift が起きないように各層の入力の分布を一定に保とうとするテクです。

具体的には、各層の入力 ( = 前層の出力 ) をバッチごとに標準化します。

 \displaystyle \mu = \frac{1}{m} \sum_{i=1}^{m} x_i

 \displaystyle \sigma^{2} = \frac{1}{m} \sum_{i=1}^m ( x_i - \mu )^{2}

 \displaystyle \hat{x} = \frac{ x - \mu }{ \sqrt{ \sigma^{2} + \varepsilon } }

ここで、 \varepsilon \sigma^{2} 0 に近い時に数値的な安定を保つためのハイパーパラメーターです。

ただ標準化するだけでは、例えば活性化関数にシグモイドを使っていた場合に線形部分しか使われなくなるかもしれないなどの問題が生じるかもしれません。

そこで、パラメータとして新しく  \beta \gamma を用意して、

 \displaystyle y = \gamma \hat{x} + \beta

と線形変換して、これを新しい入力とすれば完成です。

Batch Normalization は、Internal Covariate Shift を抑制するだけでなく、層の入力が異常に大きくなって勾配が消失してしまう問題を防いだり、学習率を大きくしたときにパラメータが発散してしまう問題を防ぐなどの効果もあるようです。

実装

gist.github.com

学習が終わって推論の段階で、バッチの組み合わせによって同じ入力が違う出力になるのがつらいときは、学習時の統計を使って  \hat{x} の変換式を固定すると良いようですが、今回は実装していません。

実験

三層のパーセプトロンで MNIST の分類を行いました。

Batch Normalization なし

f:id:joisino:20170706120414p:plain

Batch Normalization あり

f:id:joisino:20170706120428p:plain

links.BatchNormalization を使った場合

f:id:joisino:20170706120441p:plain

Batch Normalization を入れると、最終的な汎化性能は変わりませんが学習速度は速くなっています。

自作の BatchNormalization は chainer に入っている links.BatchNormalization と性能はほとんど同じになっています。

ただし、自作の BatchNormalization は自分で勾配を計算したわけではなく既存の functions を組み合わせて作ったので、速度は 4 割ほど落ちます。

最後に

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

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

参考文献

[1] Ioffe, S. et. al. (2015). Batch Normalization: Accelerating Deep Network Training by Reducing Internal Covariate Shift

Generative Adversarial Nets

Generative Adversarial Nets [1] を chainer で実装しました。

いわゆる GAN です。

最近いろいろ派生系が出ています。画像の生成モデルはほとんどこれの派生な気がします。

画像を生成する Generator (以下 G)と、画像が本物か G が生成したものかを識別する Discriminator (以下 D) という 2 つのモデルを同時に学習させていきます。

論文に書いてあるように、G は偽札製造業者、 D は警察という類推がわかりやすいです。

業者は警察にバレないように紙幣にできるだけ似せたものの作り方を学ぶ一方で、警察はより精巧な偽札を本物と区別できるように学習する、ということを繰り返すうちに、互いの技術は向上していきます。

これと同様に、G は生成するデータと教師データが D に区別できないようなパラメータを学習する一方で、 D は本物のデータか G が生成したデータかを正しく識別するようパラメータを学習していきます。

具体的には、以下の式で最適化します。

 \displaystyle min_{G} max_{D} V(D,G) = \mathbb{E}_{x \sim p_{data}(x)}[log D(x)] + \mathbb{E}_{z \sim p_{x}(z)}[log(1-D(G(z)))]

 p_x(z) は G に食わせるデータです。一様分布から生成されるノイズなどを使う場合が多いようです。

教師データを正例、G の生成したデータを負例としたときの交差エントロピー誤差をお互い最適化するということです。

論文に証明が書いてある通り、この関数の大域最適解において、G が生成するデータの分布と、教師データの分布は一致します。

あとはこれを逐次的に最適化していけばよいです。

今回の実装では、教師データに MNIST を使いました。

ソースコードは以下の通りです。

gist.github.com

生成した画像は以下の通りです。(学習の最後に出力された 10 枚を掲載します)

f:id:joisino:20170703180931p:plainf:id:joisino:20170703180935p:plainf:id:joisino:20170703180939p:plainf:id:joisino:20170703180943p:plainf:id:joisino:20170703180947p:plainf:id:joisino:20170703180950p:plainf:id:joisino:20170703180954p:plainf:id:joisino:20170703180957p:plainf:id:joisino:20170703181000p:plainf:id:joisino:20170703181006p:plain

周りにもやもやがついているのは残念ですが数字であることはわかります。

今回は G は線形 + relu + sigmoid で D は線形 + relu + batch normalization + sigmoid です。

G に batch normalization した方がよいという記事をネット上でいくつか見かけたのですが、実験してみたところ微妙でした。

以下のような感じです。

G: bn あり D: bn あり

f:id:joisino:20170703181429p:plain

G: bn なし D: bn なし

f:id:joisino:20170703181604p:plain

G: bn あり D: bn なし

f:id:joisino:20170703181537p:plain

G の bn を有効にすると数字がぼやけてしまい、 D の bn を有効にすると周辺のもやは取れるのですが学習がうまく進みませんでした。

これらの結果はハイパーパラメータの値のチューニング次第かもしれません。

今回は linear でやりましたが deep convolutional なものもまた実装したいです。

参考文献

[1] Goodfellow, I. J., et.al. (2014). Generative Adversarial Nets

Chromatic Polynomial と Acyclic Orientation

SRM717 の Hard で出題されて気になったので書き留めておきます。

Chromatic Polynomial について

まず、無向グラフの頂点から色への対応を彩色、隣接する頂点が同じ色にならないような彩色のことを正しい彩色ということにします。これはこの記事限定用語です。

以下無向グラフ  G \lambda 色で正しい彩色する場合の数を  \chi(G,\lambda) と表します

Chromatic Polynomial はこの  \chi(G,\lambda) のことです。

これは、色数  \lambda多項式になります。

まずはこの多項式を具体的に構成しましょう。

 \chi(G,\lambda) は以下の性質を持ち、また逆に、以下の性質を持つ関数は  \chi(G,\lambda) となります。

  1.  G_0 を 1 頂点からなるグラフとすると  \chi(G_0,\lambda) = \lambda
  2.  G の連結成分を  G_1,...,G_c とすると、 \chi(G,\lambda) = \Pi \chi(G_i,\lambda)
  3.  G の任意の辺  e について、  \chi(G,\lambda) = \chi(G \backslash e,\lambda) - \chi(G/e,\lambda)

ここで、 G \backslash e G から辺  e を取り除いたグラフのこととします

一方、 G / e は辺  e を縮約して、両端の頂点を 1 つに潰したグラフのこととします

以下、 \chi(G,\lambda) がこの 3 つの性質を持つことを示しましょう。

1 は、1 頂点のグラフを  \lambda 色で塗り分ける通り数なので  \lambda 通りです。

2 は、それぞれの連結成分ごとに独立に塗り分けられるので、場合の数はそれらの積になります。

3 は少し難しいです。

 e = \{u, v\} とします。

 G の正しい彩色は、 G \backslash e の正しい彩色となっています。(条件がより緩いため)

逆に  G \backslash e の正しい彩色を考えましょう。この彩色はほとんどが  G の正しい彩色となっています。

どのような彩色がダメかというと、 u v を同じ色に塗るようなもの (*) です。このような彩色はいくつあるでしょうか。

(*) を満たすような彩色  \sigma を一つ取ります。

 G / e の彩色  \sigma' として、 u v を縮約してできた頂点の色を  \sigma(u) ( = \sigma(v) ) で塗り、それ以外の頂点を  \sigma と同じ色塗りとしたものをとります。

 \sigma G \backslash e の正しい彩色であることから、 \sigma' G / e の正しい彩色であることがわかります。

また、少し考えれば、 G / e の正しい彩色  \sigma' と、 G の不正な彩色  \sigma のこの対応は 1 対 1 対応になっていることが分かります。

よって、 \chi(G \backslash e,\lambda) から  \chi(G / e,\lambda) だけ取り除いてやればよいので、 \chi(G,\lambda) = \chi(G \backslash e,\lambda) - \chi(G/e,\lambda) が成り立ちます。

以上より示されました。

逆に、3 つの性質を満たす関数があると、帰納的に  \chi(G,\lambda) と一致することが分かります。

また、以上の性質から、 \chi(G,\lambda)多項式になることも分かると思います。

Acyclic Orientation と Chromatic Polynomial について

向き付け (Orientation) とは、無向グラフの辺を向き付けて有向グラフにする操作とします。

Acyclic Orientation はその名の通り、サイクルが存在しないような向き付けのことです。

実は、Chromatic Polynomial と Acyclic Orientation は深く関係しています。

まず、 \chi(G,\lambda) に別の解釈を与えましょう。

前の章では、彩色  \sigma を、頂点から色への関数としていましたが、今回は色に番号をつけて、頂点から  \{1, 2, ..., \lambda \} への関数とします。

こうしたとき、  \chi(G,\lambda) は以下の条件を満たす、彩色  \sigma と向き付け  O の組の数となります。

  1.  O によってサイクルは生じない
  2.  O によって  u \rightarrow v と向き付けられたとすると、 \sigma(u) > \sigma(v)

正しい彩色  \sigma が与えられると、値が大きい方から小さい方に向けなければいけないので、向き付けは一意に定まります。

逆に、上記の条件を満たしているとき、辺の両端の値は異なるので、彩色は正しい彩色となっています。

よって、正しい彩色と上記の条件を満たす組は一対一対応しているので、どちらも場合の数は  \chi(G,\lambda) となります。

次に、少し変えた以下のような条件を満たす彩色  \sigma と向き付け  O の組の数  \overline{\chi}(G,\lambda) を考えます。

  1.  O によってサイクルは生じない
  2.  O によって  u \rightarrow v と向き付けられたとすると、 \sigma(u) \ge \sigma(v)

2 の条件を満たすとき  O u に適合するということにします

条件を満たす組のことを正しい組ということにします。また、彩色  \sigma が固定されているとき、サイクルを生じない  \sigma に適合する向き付けのことを正しい向き付けということにします

 G の頂点数を  N とすると、実は、 \overline{\chi}(G,\lambda) = (-1)^{N} \chi(G,-\lambda) (**) が成り立つということを示すのがこの章の山場です。

この式は、 n 個のものを m 個に分ける組み合わせの数  \binom {n}{k} と重複組み合わせ  \binom {n+k-1}{k} の関係  \binom {n+k-1}{k} = (-1)^{k} \binom {-n}{k} に相当するらしいです。

たしかに  \overline{\chi} の条件の方は重複組み合わせ感ありますが、この類推がどこまで深い対応を示しているのかは正直よく分かっていません。

それでは (**) を示していきましょう。

前の章と同様に、 \overline{\chi}(G,\lambda) を具体的に構成していきます。

  1.  G_0 を 1 頂点からなるグラフとすると  \overline{\chi}(G_0,\lambda) = \lambda
  2.  G の連結成分を  G_1,...,G_c とすると、 \overline{\chi}(G,\lambda) = \Pi \overline{\chi}(G_i,\lambda)
  3.  G の任意の辺  e について、  \overline{\chi}(G,\lambda) = \overline{\chi}(G \backslash e,\lambda) + \overline{\chi}(G/e,\lambda)

以上の式が成り立てば、示したい (**) が帰納的に成立することが少し考えれば分かります。

 G/e は頂点数が  G より 1 小さいので 3 の符号が反転していることに注意してください。

1 と 2 については前の章と全く同じです。

ここでもやはり 3 が難しいです。

 e = \{u, v\} とします。

まず、 G の正しい組について、向き付けの  e を取り除いたものは  G \backslash e において正しい組になっています。

逆を考えます。

 (\sigma,O) G \backslash e の正しい組とします。

 \sigma G の彩色にもなっていることを注意してください。

 G の向き付けとして、 O u \rightarrow v を加えたもの  O_1 O v \rightarrow u を加えたもの  O_2 を考えます。

 (\sigma, O_1) (\sigma, O_2) の少なくとも一方は  G の正しい組であって、どちらもが正しいような場合はちょうど  \overline{\chi}(G/e,\lambda) 通りであることを示します。

(i)  \sigma(u) > \sigma(v) のとき

このとき  e の向きから  O_2 \sigma に適合しません。

一方、 O_1 は適合し、サイクル  v \rightarrow w_1 \rightarrow ... \rightarrow u \rightarrow v があったとすると、 \sigma(v) \ge \sigma(w_1) \ge ... \ge \sigma(u) > \sigma(v) と矛盾が生じるので、サイクルはありません。

よって、この場合は  O_1 のみが正しい向き付けです。

(ii)  \sigma(v) > \sigma(u) のとき

この場合は (i) と同様に  O_2 のみが正しい向き付けです。

(iii)  \sigma(u) = \sigma(v) のとき

この場合は、 O_1 O_2 の両方が  \sigma に適合します。

仮に  O_1 O_2 の両方がサイクルを生じるとします。

それぞれのサイクルを  v \rightarrow w_1 \rightarrow ... \rightarrow u \rightarrow v  u \rightarrow w'_1 \rightarrow ... \rightarrow v \rightarrow u とすると、 G \backslash e O において、 v \rightarrow w_1 \rightarrow ... \rightarrow u \rightarrow w'_1 \rightarrow ... \rightarrow v というサイクルが生じてしまいます。(頂点は重複してしまうかもしれませんが特に問題ないことがわかります。)

これは、 O がサイクルを生じないという仮定と矛盾するので、 O_1 O_2 の少なくとも一方はサイクルを生じず、正しい向き付けとなります。

(iii) においてどちらも正しい向き付けとなるような場合が  \overline{\chi}(G/e,\lambda) 通りであることを示します。

 (\sigma,O_1) (\sigma,O_2) がどちらも  G において正しい組となるような  (\sigma,O) G / e の正しい組は一対一に対応することが以下のように分かります。

 G / e の彩色  \sigma' として、 u v を縮約してできた頂点の値を  \sigma(u) ( = \sigma(v) ) とし、それ以外の頂点を  \sigma と同じ値としたものをとります (***)。

 e をどちらの向きにしても  O G でサイクルは生じないので、  O G / e でサイクルを生じません。

また、 O G \backslash e \sigma に適合してるので、 G / e では  \sigma' に適合しています。

よって、 (\sigma', O) G / e において正しい組であることがわかります。

逆に、 G / e の正しい組  (\sigma',O) をとり、 \sigma' に(***) の逆の対応をさせた  \sigma を考えると、 O G \backslash e \sigma に適合しています。

 O G / e でサイクルを生じないので、 e をどちらの向きにしても  G でサイクルは生じません。

よって、重複して数えなければならない場合が  \overline{\chi}(G/e,\lambda) の要素と一対一に対応し、 \overline{\chi}(G,\lambda) = \overline{\chi}(G \backslash e,\lambda) + \overline{\chi}(G/e,\lambda) となることが分かります。

以上より示されました。

 \lambda = 1 のときを考えると、向き付けが正しいことと、サイクルが生じないことが同値なので、無向グラフの Acyclic Orientation の数は  \overline{\chi}(G,1) = (-1)^{N} \chi(G,-1) となることが分かります。

SRM 717 Div1 Hard AcyclicOrientation

無向グラフが与えられるので Acyclic Orientation の数を mod 6 で求めてくださいという問題です。

中国剰余定理より、 mod 2 と mod 3 で  \overline{\chi}(G,1) = (-1)^{N} \chi(G,-1) が求まればよいです。

また、 \chi(G,\lambda)多項式であることを考えると、引数も mod を取ってしまってよいことが分かります。

よって、求めるものは mod 2 で  \chi(G,1) と mod3 で  2^{N} \chi(G,2) です。

前者は辺があれば 0, 辺がなければ 1 です。

後者は、 \chi(G,2) は 2 彩色の場合の数なので、二部グラフでなければ 0, 二部グラフであれば  2^{N} \chi(G,2) = 2^{N+C} ただし  C は連結成分の個数です。

あとはこれらをループを回して 6 通り調べるなどして mod 6 に戻してやれば解けたことになります。

最後に

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

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

参考文献

[1] Stanley, R. P. (1972). ACYCLIC ORIENTATION OF GRAPHS.

計算機科学実験及演習3ハードウェア(CPU製作)記

先日、長いようで短かった実験がついに終わったので記録を残しておきます。

はじめに

京都大学情報学科の計算機科学コースでは計算機科学実験及演習(以下、実験)という必修科目が 2 回生前期 / 後期、 3 回生前期 / 後期の計 4 つあります。

実験 1 はプログラミングの入門みたいな内容
実験 2 はマリオ AI の作成と電子回路
実験 3 は CPU 製作とインタプリタ製作
実験 4 はいくつかの分野から選択

という感じです。

僕たち 2015 年入学の学生から実験の内容が少し変わったのですが、実験 3 の CPU 製作は昔からあるみたいです。

CPU を製作する学生実験といえば某大学某学科の CPU 実験が有名で知っている人も多いと思います。

CPU 実験でググるとすごい人の製作記がたくさん出てきて面白いのですが、いまググるとこの記事が霞んで見えるので読み終わってから調べてください(参考文献にいくつか並べてあります)

あちらの実験は有名なのに、実験 3 の記事は(調べ方が悪いのか)調べてもあまり出てきません。

これではあまりにも寂しいし、僕自身実験をはじめる前に製作記みたいなのを読みたかったので後輩のためにも少しがんばって書くことにしました。

実験について

実験 3 (ハードウェア)では二人がかりで FPGA を使って CPU を製作します。

期限は二ヶ月弱。今年は初回授業が 4/13 で、デモが 6/2 でした。

実験の時間は毎週木曜 3 限から 5 限と金曜 1 限から 5 限までの 8 コマ(毎週 12 時間)です。めちゃくちゃ長い。

しかも 4 分の 1 以上休むと単位がもらえないというスパルタぶりです。

作るアーキテクチャについては SIMPLE という 16 ビットアーキテクチャが提示されて、これに準拠したものを作ります。

準拠するといっても改良するのはアリで、どこまで変えていいのか聞いてみたところ面白ければなんでもいいらしいのでほぼ自由です。

ただし、レジスタ幅/データバス幅は 16 ビットでないとダメです。この制約が割と辛いので注意してください。32 ビットアーキテクチャにするくらいなら許してほしいところです。

この制約はソートコンテストで SIMD の幅増やす勝負とかにならないようにするためだそうです(それはそれで面白そうですが)

ソートコンテストというのは、製作した CPU 上で 1024 個の 16bit の符号付き整数をソートする速度を競うコンテストのことです。

必須参加ではありませんが、だいたいこのコンテストに参加して上位を狙うのを目標にする班が多いようです。


とりあえず一通り説明したので、これからは日記形式でやったことを説明していきます。

日記

4/13

実験初日です。

ハードウェアにはめちゃくちゃ苦手意識があるので気が重かったです。

初週はまだ CPU を作らないようで、HDL の説明を軽くされたあと、導入課題として個人で 10 進カウンタを作って 7 セグ LED に表示するというものが出ました。

カウンタ自体はそう難しくないんですが、 LED への出力の仕様がよく分からずハマりました。

こういうハードウェア部分でハマって、マニュアルを見ても分からないときは早めに TA に相談するのが良いと思います。

▲ハードウェア難しすぎて絶望する joisino 氏

4/14

まだカウンタが終わらないので作り続けます。


昼になった時点で CPU 製作の班分けが発表されました。

僕のペアは kmcyuma くんでした。

KMC の競技プログラミング勢ということで、気心が知れて、プログラミングがとてもできる人なので助かりました。

ペアはたいへん重要なのですが、自分たちで選べる訳ではなく、指定されます。

これから実験 3 に挑む人は良い人と当たるために春休みのうちに徳を積んでおきましょう。


午後は、SIMPLE アーキテクチャの説明と git の使い方説明会がありました。

今年は全員 HDL(Verilog HDL か VHDL) で CPU を作ることになったのですが、去年までは HDL で作る班と CAD でポチポチ配線する班があったそうです。

CAD で CPU 配線するの激ヤバだと思うんですが、今年は全員 HDL ということで良かったです。

僕は Verilog HDL よりも System Verilog の方が良かったのでそれでもいいかと聞いたら OK が来たので System Verilog で作ることにしました。


なんとかこの日のうちにカウンタを作り終えて、レポートを書き終わるところまでたどり着きました。

レポート部分が意外と重いので注意してください。


最後の 1 時間くらいは相方と今後の CPU 製作について打ち合わせしました。

だいたい次のようなことを相談しました。

  • 相方はシミュレータとアセンブラとソートプログラムを作る。
  • 僕は CPU を作る。
  • とりあえず来週くらいには単位コアを作って精神を安定させる。
  • 単位コアはすぐに捨てていい感じのコアをつくりはじめる。
  • さすがにパイプライン化はしたい。
  • OoO とかスーパースカラ化もできればやりたい。
  • ソートコンテストで 1 位になるのを目標にする。
4/18

CPU が完成しました

製作パートどこいったんだよ製作記じゃなかったのかよという感じですが僕自身もビックリするくらいすぐ完成してしまいました。

マイクロアーキテクチャの図を見ながら、ひとつずつ部品を作って配線していくと出来上がります(それはそう)

コアを作る上でオススメなのは

  • SIMPLE/B というマイクロアーキテクチャが参考に上がっていますが、あれよりもパタヘネに載ってるやつのほうが慣れてて書きやすいです(これは人によりそう)
  • コンパイルに 2 分くらいかかってつらいですが、Modelsim での RTL レベルのシミュレーションなら 1 秒くらいでコンパイル終わるのでこちらでバリバリデバッグして、最後にタイミングの確認とかの段階になって論理合成するのが良いです。
  • 機能を増やすにつれてどこに何があるのか混乱してくるので、コアの全体図を描いておくとわかりやすいです。
  • テスト用にいろいろな難易度のプログラムを書いておくと良いです。僕たちの班は入力をそのまま出力に返す、演算して出力する、1 から 10 までの和をループで求める、フィボナッチ数列の 10 番目を再帰で求める、バブルソートをするなどをテスト用に書いて、ハードがバグったら簡単なものからどこまで通るかを確認してバグ箇所を特定したりしていました。

特にふたつ目は生産性爆上がりなので、ぜひ実行してみてください。

ちなみにアセンブラとソートプログラムも相方がまだ作っていなかったので僕が作ってしまいました。

愚直なクイックソートで 4.96ms でした。

▲完動を喜ぶ joisino 氏

4/20

5 段パイプライン化しました。

パタヘネに有名な文言「誤信: パイプライン化は容易である」というのがありますが、意外とすんなりパイプライン化できました。

ここで特にパタヘネのマイクロアーキテクチャを参考にしたのが効いてきたのだと思います。

SIMPLE/B ベースでパイプライン化しようとするとハザード検出を自分で考えることになるのでバグらせやすいですが、パタヘネのマイクロアーキテクチャは本に書いてある上に講義でもパイプライン化の手順をやったはずなので、手際よくできるはずです。

フォワーディングを実装して、この日にソートコンテストに一度目の提出をしました。1.451ms でした。

4/21

そろそろ単位コアを捨てて新しいコアを作り始めてもよかったのですが、相方と話し合った結果、このままでも歴代最速を更新できるんじゃないかという気がしてきたのでもう少し作り続けようという結論になりました。

ちなみに歴代最速は 2012 年度の 4 コアプロセッサで 0.3680ms です。

4/23

分岐予測を実装、クロックを上げる、Jump And Link 命令すらなかったので追加するなどなどを行って歴代最速を更新しました

小さい時は挿入ソートにする & ピボット選択をがんばるクイックソートで 0.3650ms です。(これは提出していません)

もうこれ以上の高速化は限界っぽかったのでマルチコアにしようかということになってきました。

そういえば、このあたりは他の講義中もずっと CPU 作っていて、そのせいか週 3 くらいで CPU を作る夢を見ていた記憶があります。

特に CPU の配線を全て手作業で行わされる悪夢にうなされたのが印象に残っています。

4/27

基数ソート + 挿入ソートを書いて 0.1916 ms で提出しました。

アルゴリズムの変更でここまで速くなるのはアルゴリズムの勝利という感じがして嬉しかったですね。

twitter で自慢する joisino 氏

4/28

8 コアプロセッサが動いて 0.05040ms で提出しました。

▲記録の更新に浮かれる joisino 氏

結果的にこれがソートコンテストの最終提出になりました。

正直 4 月中にマルチコアが動くと思ってなかったので動いた時はめっちゃ嬉しかったです。

俺は 2 週間でマルチコアプロセッサを作ったぞってあと 10 年くらいは自慢するかもしれませんが聞き流してもらって構いません。

ちなみにメモリとかタイミング制約を限界まで取って 8 コアでコンパイルすると論理合成に 1 時間くらいかかって絶望するのでそこそこで我慢しましょう。6 割にするだけで 5 分くらいになります。

マルチコアですが、ソートコンテストが動く最低限の同期だけ実装したものならそこまでヤバくはありません。

僕の班でも、8 つのコアはほとんど独立で、メモリだけ共有しています。

動作するために特別やったことといえば、メモリの排他制御が必要なのでセマフォ入れたぐらいです。

あとは高速化のためにメモリアクセス権のプライオリティをクロックごとにコアに順番に割り当てていくことと、メモリの 2 ポート化(読み込み 1 / 書き込み 1)をしました。

ソートコンテストでは要素数 1024 個の配列をソートするのですが、ハードウェアの制約上メモリは 2 ポートまでしか実現できず、クロック数も 250 MHz が上限(といっても CPU を 250 MHz で動かすのは至難の業ですが)なので、仮に値の収納先がわかって毎クロック値を読み込み / 書き込みを行ったとしても 1024 / 250 MHz = 0.0041 ms くらいはかかります。

相当無理をした仮定でも 13 分の 1 くらいにしかならないことを考えるとこれ以上の本質的な高速化は厳しいという結論になり、これからはソートコンテスト以外の道を模索することになりました。(後から思うと、僕たちの CPU は 90MHz なので、パイプライン段数を増やしてクロックを上げればまだ高速化はできたと思います。)


ところで少し話は変わりますが、この実験では中間レポートというものがあり、我が班は非常に焦っていました。

というも、実験の評価項目では分担がうまくできているかが重視されているのですが、ここまで CPU は全て僕が作ってしまっていました。

ということで、この日に急遽相方に CPU の仕組みと System Verilog を教えて、細かい部品を作ってもらうことにしました。

4/29

ゴールデンウィークですが、コンパイラを作れたらなぁと思って LLVM のドキュメントを眺めたりしつつ遊んでいました。

あとゴールデンウィーク中にちょくちょくレポートを書きましたがかなり分量が多いので気をつけてください。

5/6

LLVM バックエンド難しすぎて挫折する。

5/7

とりあえず flex + bison で温かみのある昔ながらのコンパイラを作ることにしました。

ソースは C サブセットにするつもりだったのですが厳密にサブセットにするのは面倒だったので C サブセットっぽい何かになってしまいました。

この日は演算命令をコンパイルするくらいを実装しました。

5/8

コンパイラに比較、if、else、whileあたりを実装しました。

このあたりで気づいたのですが、16 ビットということもありアーキテクチャが弱すぎでコンパイラを作るのが相当大変です。

たとえば、分岐命令は上下 128 命令くらいまでしか飛べないので if 文などでそれ以上飛ぼうと思うと一旦別の場所に飛んでから遠くの場所に飛ぶ必要があります。

ここまではよくあるテクなのですが SIMPLE は即値分岐や即値代入も絶対値 128 くらいまでしか指定できないので困り果てました。

調べてみると MIPS は即値分岐で飛べない場所に行かないように配置するらしいのですが流石に 128 しか飛べないように配置するのでは制限が大きすぎます。

ラベルの x ビット目から y ビット目を表すラベルみたいなのを作って無理やり処理したのですがこういうのはどのように処理するのが正解だったのでしょうか。


あと、この日にようやく相方が作っていたシミュレータがまともに動くようになりました

もう既に完動してしまっていますが、いちいち波形シミュレータを使うのは面倒すぎるので、このシミュレータはこれからのコンパイラ製作に大いに役立ちました。

5/9

コンパイラに関数呼び出しを実装しました。

5/10

コンパイラに関数の引数や配列を実装しました。

このあたりで以下のフィボナッチ数を求めるプログラム程度ならコンパイルできるようになりました。

gist.github.com

5/11

中間レポートを提出しました。

相方もいい感じにパーツを作ってそれらを CPU に組み込んでなんとか協力した感じを醸し出しました。

5/12

この日は中間デモでした。

完成品を見せるというよりは進捗確認の意味合いが強い感じで、とりあえず 1 命令以上実機で動くことを見せるというものでした。

張り切った僕たちの班はすでにコンパイラの原型までできていたので、ちょっと C プログラムを書いてマンデルブロー集合を出力するデモをしました。

f:id:joisino:20170604175836p:plain
▲誤差がすごいマンデルブロー集合

FPU とかはないので固定小数の掛け算のライブラリを書いて動かしました。


とりあえず山場の中間デモまで終わったので最終デモでやる余興を何にするかという話になりました。

いろいろ考えてはみたのですがハードウェアの制約がかなり厳しいです。

FPGA からなにか使える端子が伸びていてディスプレイでも繋げられたらいろいろやりようがあると思うのですが、残念ながら出力は LED とブザーしかありません。

OS っぽいものを作る案もあったのですが、ハーバードアーキテクチャで作ってしまったのでプログラムの展開をどうするんだという話になって流れてしまいました(後から思うとハーバードでも何かしらは作れたと思うのですが)

とりあえずボードに乗っているハードウェアを最大限使ってテトリスを作ることにしました。

テトリスのプログラムは相方が作って、ハードウェアとコンパイラ方面を僕が作ることにしました。

とりあえず僕は LED の出力周りをいい感じにアーキテクチャから見えるようにしたりしました。

f:id:joisino:20170604175957j:plain
▲実験で使うボード

5/15

コンパイラに for や continue や break などを追加したりしていました。

5/19

ブザーで音楽を鳴らせるようにしました。

音楽部分はプロセッサから独立させてひたすらメモリからデータを読みだして BGM を鳴らし続けるようになっています。

アーキテクチャから見えるようにしたかったのですが、BGM を鳴らすとなるとタイマ割り込みとかが必要になってきそうなので諦めました(これも後から思うと色々やりようはあったと思います)

だいたいテトリスが動くようになりました。

割と大変だったのが、メモリの容量不足です。

実験で使うボードでは 16 キロワードくらいまでしかメモリが入らなくて、データ 8k 命令 8k にしているのですが、コンパイラが最適化をしてないこともあり、普通にテトリスを書くと溢れてしまいました。

なんとか相方に手で命令数を最適化したりしてもらっていました。

5/25

この時点になってもまだソートコンテストの参加者が僕たちだけだったので、まだ完成していない KMC 部員の席に行っては煽ったりしてました(ごめんね)

5/26

タイマーを追加して、リセットからの時間を取得できるようにしました。

他にもいろいろと微調整をしてデモ用のプログラムが完成しました。


さて、このあたりでまたまたレポート問題です。

最終レポートはこれまでのレポートよりも分量が多く、班で提出するレポートが 4 種類と個人で提出するのが 1 種類あります。

なかでもユーザーズマニュアルというのが大変で、「IPコア(ソフトコアプロセッサ)として提供することを想定し、プロセッサを使用するうえで必要十分な情報を記載すること」といういかにも大変そうな注意書きが書かれています。

特に僕たちの班はごちゃごちゃと機能を付け足してしまったので書くことが山盛りです。

そんなわけでこのあたりから実験中はずっとリファクタリングをしながらレポートを書き続ける回になりました。

6/1

デモ前日です。

デモはペアがそれぞれ絶対に喋らないといけないのに、相方はハードウェアの中身については知らないらしいのでいろいろとレクチャーして発表の分担を考えたりしました。

6/2

いよいよデモ当日です。

デモといってもスライドを作ってみんなの前で発表するわけではなく、TA と教員が班席まで回ってきて 10 分くらいプロセッサやプログラムの説明をする感じです。

知り合いの発表は聞きにいったりしていましたが。

テトリスはこんな感じに仕上がりました。

youtu.be

声が入っているので音は入れていませんが、プレイ中はコロブチカが流れます。

あと相方がプレイして僕が撮ってるのですがミスするたびに笑ってカメラを揺らしています。ごめんなさい。

これだけ色々言ったわりには大したことないようなという感じかもしれませんが、とにかくデモが終わった後はようやく終わったという開放感でいっぱいでした。

レポートもこの日のうちになんとか書き終えて提出することができました。

6/4

この記事を書いています。

感想

2 ヶ月弱、つらいことも多かったですが、振り返ってみると全体的に楽しかったと思います。もうやりたくはありませんが。

期間も短く班員も 2 人なのであまり大層なことはできませんが、こんな短期間で CPU 製作を一通り学べるのはものすごく良い経験だと思います。

僕たちの班も最後の方は若干だれてしまったので、全力疾走するにはこの位の期間が限度なのかもしれません。

色々やりたいことの中から、無理だと諦めたものもありましたが、最終的に作れたものには満足しています。

来年以降実験を受ける人たちは僕たちが諦めたものたちも含めて、僕たちの成果物なんかよりもっとすごいものを作ってください。(半分そのために記事を書いているので期待してますよ!)

的を射ているかは分かりませんが若干アドバイスです

  • OoO とスーパースカラ

- ソートコンテストをやるうちにいつの間にか目標から外れていました。というのは、通常実験ではアセンブリコードを手書きするので、相当のハザードを手で解消できてしまい OoO の恩恵をあまり受けられないというのと、そもそも 8 コア化するとすごい量のメモリアクセスがあってアクセスキュー待ちがボトルネックになるのでこのあたりの高速化テクの意味が無くなってしまったからです。ですが、正直このあたりの技術を実装するのを目標にしている人たちならソートコンテストは無視してもっと大規模なプログラムを目標に据えて、OoO とかスーパースカラを実装していくのがいいと思います。

  • second architecture

- 僕たちの班も最初は単位コアを作って新しいアーキテクチャに移ろうと考えていたのですが、結局最後まで first を建て増しし続けて終わりました。これについてはどちらが良いとも言えないと思います。もうちょっと実験期間が長ければ有用だと思うのですがなにしろ実験期間が短いので second を作る余裕はあまりありません。second を作ることにするならちゃんとペースを考えたほうが良いと思います。

  • OS

- OS を作りたいなら最初の設計の段階から決めてかかったほうが良いです。ハーバードで作るよりノイマン型で作ったほうがいろいろやりやすいと思うし、割り込みとかもあんまりすぐに建て増しできる感じではないので最初の設計で考えたほうがよいはずです。たぶんハードウェアに詳しくない人は最初からこのへん考えるのは難しい(少なくとも僕はそうだった)ので、OS 作ったりする場合は second 作るのが有用なのかもしれません。

  • 出力

- ボードに出力が全然ないと言いましたが、唯一出てるリッチな出力がコンフィギュレーション用の USB Blaster です。全く知りませんが、これを使っていい感じにシリアル通信したりできるかもしれません。最後の方になって気づいたのですが結局調べる時間もないまま終わりました。これが使えるならできることが一気に広がると思うので調べてみると良いかもしれません。


レポジトリを置いておきます(剽窃があるといけないので公開しないようにと先生が言っていましたが、頼んだら許可してくれました。来年以降受講する人は絶対に剽窃しないでくださいね。)

github.com

一息つく間もなく来週からはソフトウェア実験がはじまるのでがんばります。

参考文献

ためになった本たちを紹介しておきます。

CPUの創りかた

CPUの創りかた

めっちゃいい本。あと絵がかわいい。

今回の実験で使ったわけではなく昔読んだものですが、この本のおかげで CPU がどんな風にできているのか理解できた記憶があります。

フリップフロップのイメージはまだこの本に書いてあるカメラのページを思い出します。

パタヘネとか読んだことあるならもう読む必要があるかは分かりませんが、ハードウェアがめちゃくちゃ苦手な意識のある人は読んでみるといいかもしれません。

コンピュータの構成と設計 第5版 上

コンピュータの構成と設計 第5版 上

コンピュータの構成と設計 第5版 下

コンピュータの構成と設計 第5版 下

いわゆるパタヘネ。

講義とかで使ったことがあるはずです。相方はなぜか持ってませんでした。

これがないとさすがに実験 3 を乗り越えるのは難しいのでちゃんと買って読みましょう。実験 3 では特に上巻が重要です。


ディジタル回路設計とコンピュータアーキテクチャ[ARM版]

ディジタル回路設計とコンピュータアーキテクチャ[ARM版]

ディジタル回路設計とコンピュータアーキテクチャ (IT Architects' Archiveクラシックモダン・コンピューティング)

ディジタル回路設計とコンピュータアーキテクチャ (IT Architects' Archiveクラシックモダン・コンピューティング)

パタヘネに比べるとあまり有名ではないかもしれません。

僕は上の ARM 版を読みました。

概念自体はパタヘネを読んだことがあるならあまり新しいことは無いですが、この本の良いところは System VerilogVHDL (通常版は Verilog HDL と VHDL らしい?)のサンプルコードが載ってあることです。

最後にはシンプルなシングルサイクルの ARM プロセッサのコードまで載っています。

僕は System Verilog をだいたいこの本で学びました。

HDL を学べる資料は貴重なので、実験初期に HDL がよく分からなかったら読んでみると良いかもしれません。

コンピュータ設計の基礎

コンピュータ設計の基礎

高性能コンピュータ技術の基礎

高性能コンピュータ技術の基礎

kindle でセールしてたので買いました。

上のコンピュータ設計の基礎については、パタヘネに書いてあることがほとんどです。

下の高性能コンピュータ技術の基礎は、OoO とかスーパースカラがかなり丁寧に解説してあります。

OoO とかスーパースカラを実装しようかと思っていて、ヘネパタは高い/重いと思う人は読んでみるとよいと思います。

ヘネシー&パターソン コンピュータアーキテクチャ 定量的アプローチ 第5版

ヘネシー&パターソン コンピュータアーキテクチャ 定量的アプローチ 第5版

いわゆるヘネパタ。

OoO とかスーパースカラあたりの部分だけ読みました。

実験で使うとしたらそのあたりだと思います。

僕もちゃんと読めてないのでまた時間があるときに挑戦したいです。


nullpo-head.hateblo.jp
kw-udon.hatenablog.com
sueki743.hatenablog.jp
github.com

某大学某学科の CPU 実験のすごい人たちの記事(一部)です。

ほかにも「CPU実験」とかで調べるとたくさん出てきます。

この実験がはじまる前〜初期くらいはいろんな製作記を読みまくっていました。

読み物としても楽しいし、モチベも上がるのでオススメです。



とても長い記事になってしまいましたが、最後まで読んでいただきありがとうございました。

最小全域有向木(m log n)

明日、KMC の組合せ最適化読書会でやる予定なのですが、少し調べた限りでは日本語の解説が見つからなかったので書いてみます。(注:組合せ最適化本体には O(nm) のアルゴリズムしか載っていません。)

 

こういうアルゴリズムに慣れていない人向けに書いたつもりなのでプロ各位にはまどろっこくて分かりづらい表現になったかもしれません。ゆるして。

最小全域有向木について

辺にコストのついた有向グラフが与えられて一つの頂点が指定されたとき、その頂点を根とする全域有向木でコスト最小のものを求める問題です。

 

根用の頂点を追加して、各点にコスト INF の辺を張れば頂点を指定しない最小全域有向木が解けます。コスト 0 で張ると最小有向森です。

 

有名な解法に Edmonds / Chu-Liu のアルゴリズムがあります。この計算量は O( nm ) で、組合せ最適化や Spaghetti Source - 最小全域有向木(Chu-Liu/Edmonds) に載っています。今から説明する O( m log n ) の解法は Tarjan [1977] によるもので、このアルゴリズムの亜種(というか高速化しただけ)です。

アルゴリズム

  1. 根以外の頂点を適当に選びます
  2. 今いる頂点 v に入ってくる(自己ループでない)辺のうちコストが一番小さい辺 ( u -> v ) を取ります。今いる頂点が圧縮された頂点なら、もとのサイクルで v に入ってくる辺を使う辺から除きます。u に移り、( u -> v ) を使う辺に追加します。
  3. 使う辺でサイクルができたら、サイクルを一点に圧縮します。このとき、サイクル外からサイクル上の頂点 w に入ってくる辺のコストを、w に入ってくるサイクル辺のコストだけ減らします。
  4. まだ根と繋がらなかったら 2 に戻ります。繋がったときは、他に頂点が無いなら終了。まだあるなら適当に頂点を選んで 2 に戻ります。

やることはすごく簡単です。少し下で解説をしているのでこれを読んで分からなければそっちを読んでみてください。(むしろこれだけで分かったらすごい)

解説

有向木では、根以外の入次数は 1 です。(もちろん根は 0 です。)なので、各頂点について入ってくる辺をちょうど一つ選ぶわけです。このアルゴリズムでは、手順 2 で貪欲に一番コストが小さい辺を選んでいっています。サイクルができなければこれが一番良いことはすぐに分かると思います。

 

厄介なのはサイクルができたときです。サイクルができた時は、そこで詰んでしまうので、サイクルを「開いて」やらなければいけませんf:id:joisino:20170111142108p:plain

サイクルに繋がる辺を選ぶ訳ですが、上の図のように開くと総コストが d - c だけ増えます。(コスト c の辺を使うのをやめてコスト d の辺を使うため。)このアルゴリズムではこのコストについて貪欲に小さいものを選びます。アルゴリズムの手順 3 で「このとき、サイクル外からサイクル上の頂点 w に入ってくる辺のコストを、w に入ってくるサイクル辺のコストだけ減らします。」と言っているのはそういうことです。(「サイクル外からサイクル上の頂点 w に入ってくる辺のコスト」は図の d に対応し、「w に入ってくるサイクル辺のコスト」は図の c に対応します。)

 

サイクルができるたびに上のように対応して、サイクル(だったもの)はそれ以降変更しないように一点に圧縮してしまいます。(これは、正当性の証明のところで示すようにサイクルの一辺以外を全て含むような最小全域有向木があるので OK です。)

 

入ってくる辺を貪欲に選んで、サイクルを貪欲に切って圧縮するというのは Edmonds のアルゴリズムも同じです。違うところは、Edmonds では全部の頂点について入ってくる辺を選んでから圧縮するのを繰り返しているのに対し、このアルゴリズムでは DFS の要領で一つずつ頂点を処理していくことです。

 

現在の頂点から辺を逆に辿ってパスが伸びるか、サイクルができても圧縮してしまうので、現在の処理している辺たちは常にパスのように見えることになります。

 

こうすると、辺を辿るたびに、パスが 1 つ伸びるか、サイクルができて、サイクルの大きさくらいの頂点が消えるので、 2n ステップくらいで止まることが分かります。

正当性の証明

Edmonds のアルゴリズムの正当性の証明とほとんど同じです。そっちを知っている人はこれも証明できると思うので自力でやってみてもいいかもしれません。

 

証明は圧縮する回数についての帰納法で行います。

 

圧縮する回数が 0 回(圧縮しない)とき

これは、解説のところでも言及しましたが、ほとんど明らかだと思います。

このアルゴリズムで得られた有向木 T と、任意の有向木 S を比べたとき、各頂点に入ってくる辺同士についてコストを比べると、T の辺は入ってくる辺の中で一番コストが小さいものを選んでいるので、もちろん T の辺のコストが S の辺のコスト以下になります。ということは、総和をとっても T 全体のコストは S 全体のコスト以下です。

 

圧縮する回数が k + 1 回のとき

まず、最初に圧縮することになるサイクルを C とします。 C は L 本の辺からなるとします。

 

ここで、なんとサイクル C の辺のうち、L - 1 本の辺を使う最小全域有向木が存在することがいえます。(この証明は組合せ最適化にも載っています。慣れない人には少し難しいかもしれないので分からなかったら認めてしまって次に進んでも良いかもしれません。飛ばす場合は「補題の証明終わり」というところまで進んでください。)

 

背理法により補題を証明します。

最小全域有向木のうち、サイクル C の辺を使う本数が最も多いものを A とします。C の辺のうち K ( >= 2 ) 本使わないがあると仮定します。

C の辺のうち、A で使われていない辺を、一周回る順番に ( a_1, b_1 ), ..., ( a_K, b_K ) とします。ここで、i = 1, 2, ..., K について A には b_i から b_{i-1} にパスがあることを示します。( i = 1 のときには i-1 は K を表していると思ってください。)

( a_i, b_i ) が使われていないのだから、A には b_i  に入ってくる他の辺 e があるはずです。

e を使うのやめて ( a_i, b_i ) を使うことにしてみます。

( a_i, b_i ) のコストは b_i に入ってくる辺の中で一番小さいのでコストの総和は変更前以下になります。変更後もちゃんと有向木だったとすると、コストが真に小さくなっても A が最小であることに反するし、同じでも C のうち使う辺の本数が 1 本多いので、 A のとり方に違反します。

考えられるのは使う辺を変えたことでサイクルができてしまったということです。

これは、A において、b_i から a_i にパスがあったということです。a_i からそのパスを逆向きにたどって行くと以下の図から分かるように b_{i-1} にたどり着きます。(入次数が 1 なので一意にたどれます。また、(a_{i-1}, b_{i-1}) は使わない C の辺のうち一つ前の辺なので、間の辺は全て A でも使われています。)

f:id:joisino:20170111211948p:plain

つまり、A 上には b_i から b_{i-1} にパスがあるということです。(b_i から a_{i-1} へのパスの一部分。)

これでパスの存在は示せました。

しかし、これは b_K からb_{K-1} にパスがあり、b_{K-1} から b_{K-2} にパスがあり、...、b_1 から b_K にパスがあるということになり、A にサイクルがあるということになります。

これは A が有向木であるという仮定に反します。

よって、最小全域有向木のうち、サイクル C の辺を使う本数が最も多いものは L - 1 本の辺を使うということになります。

補題の証明終わり

 

もとの証明に戻ります。圧縮する回数が k+1 回のときにアルゴリズムで得られた有向木が最小全域有向木であることを示すのが目的でした。

 

このアルゴリズムで得られた有向木を T とし、サイクル C を圧縮したときにできるグラフ(以下圧縮グラフと呼びます)で T に相当する有向木を T' とします。帰納法の仮定より、 T' は圧縮グラフでの最小全域有向木です。

補題により存在が示された、C の辺のうち L - 1 本を使う最小有向木を T* とし、圧縮グラフで T* に相当する有向木を T*' とします。

ここで、 T のコストは T' のコストにサイクル C に含まれる全ての辺のコストの和を足したものとなります。T* についても同様に T*' のコストに C のコストを足したものになります。

これは、 T' で C に相当する部分以外のコストは同じで、 C は、1 本以外全て使うことになりますが、使わない辺については C に入ってくる辺のコストの変更によって調整されているので、このようになります。

よって、( T のコスト ) = ( T' のコスト ) + ( C のコスト ) <= ( T*' のコスト ) + ( C のコスト ) = ( T* のコスト )

となり、 T が最小全域有向木であることが分かります。真ん中の不等号は T' が圧縮グラフ上での最小全域有向木であることからです。

f:id:joisino:20170111225913p:plain

 以上、帰納法によって正当性が示されました。

 

計算量

計算量を調べます。

解説のところで言及したように、1 ステップ ごとにパスが伸びるか、頂点がマージされて減るのでアルゴリズムの手順 2, 3 を実行する回数は O(n) 回です。

 

各頂点について入ってくる辺を隣接行列でもって愚直に管理するとします。手順 2 で最小値を求めるのは毎回 O(n) で、これは O(n) 回行われます。手順 3 で圧縮するのも、1 つの点につき O(n) でマージできて、合計高々 n 頂点しかマージされないので、結局全体の計算量は O(n^2) になります。

密グラフでこのアルゴリズムを走らせたいときはこれ以上速くなりようがないのでこれで終了です。

 

次に疎グラフで O( m log n ) を達成するやり方を考えます。

dijkstra 法などを学んだ人ならなんとなく分かると思いますが、辺をヒープ(prioirty queue)で管理します。

こうすると手順 2 で最小値を求めるのは毎回 O(1) になります。

 

問題は圧縮です。ヒープをうまくマージしたりコストを一様に減らしたりする方法なんてあるのかと思いますがあります。

一つには、マージテクを使うとできます。これで O( m log n + n log^2 n ) とかになります。

もう一つには、meldable heap というものを使います。これはマージが効率よくできるヒープのことで、leftist heap や skew heap などがあります。blog - Leftist Heap & Skew Heap に詳しいのでそちらを参照してみてください。これで O( log n ) でマージできます。また、コストを一様に減らす方法ですが、ヒープのノードに一様に足された数を持っておいてマージするときに遅延伝播すればうまくいきます。(segtree の遅延伝播と全く同じです。)これで全体で O( m log n ) になります。おめでとうございます。

 

具体的な実装例は使い道のところを参照してみてください。

 

ちなみに、fibonacci heap を使うとマージが O(1) でできます。コストを一様に減らす解析は読んでないので良く分かっていませんが Gabow et al. [1986] によると全体で O( m + n log n ) になるらしいです。誰か読んで分かりやすく教えて。

使い道

競技プログラミングでは今のところあまり実用性はないかもしれませんが、とりあえず首都(The Capital | Aizu Online Judge)をアルゴリズムパンチで殴り倒せます。

 

実装例です。多めにコメントを書きました。MSA という関数が、 r を根とする最小全域有向木のコストを求める関数です。「ここが本質です。」と書いてあるループ部分が最小全域有向木のメインの部分です。

The Capital ( AOJ 2647 )

 

終わりに

この記事を書くまでは最小全域有向木は無向の最小全域木よりも(少なくとも競技プログラミングで使うようなアルゴリズムでは)有意に計算量が大きいものだと思っていたんですが、割と簡単なアルゴリズムで似たような計算量が達成できるのにビックリしました。

最小全域有向木を使う問題はあんまり見ないですが、n = 10^5 とかで割と簡単に書けるので使いようがあるんじゃないかと思いました。(浸透しないとただの知識ゲーだけど)

 

 

ここが分からないとかここ間違ってるよとかあれば気軽に言ってください。

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

 

(1/11 23:33 マージテクでも少しいい感じにできるよということなど少し追記しました)

参考文献

[1] B.コルテ, J.フィーゲン [2012]: 組合せ最適化 第二版. 丸善, 2012, 第 6 章.

[2] Gabow, H.N., Galil, Z., Spencer, T., and Tarjan, R.E. [1986]: Efficient algorithm for finding minimum spanning trees in undirected and directed graphs. Combinatorica 6 (1986), 109-122

[3] Tarjan, R.E. [1977]: Fining optimum branchings, Networks (1977), 25-35

2015ふりかえりと2016めも


昨日振り返りするの忘れていたので元日に振り返りと目標どっちも書こう

 

2015年の振り返り
こんなことを買いていたのでその振り返り

 

大学に入る
はい

 

赤くなる
TC 1584 -> 1827 (+243)
CF 1876 -> 2150 (+274)
うーんダメ
一瞬だけCF赤くなったけど色基準変更されてレートも落ちてしまった
来年はちゃんと精進して赤くなるぞ

 

ゆゆ式二期が決まる
ダメですね......

 

良い一年だったけど目標は全然達成できてないですね
2015年は全然競プロできなかったので今年は精進するぞ

 

2016年めも
TC,CF赤くなる
ゆゆ式二期が決まる