Notebook

Flask+Nginx+uWSGI+CentOS7でアプリケーションを作る際にやったこと+日記::12/9

ConoHa1のCentOS7上でFlaskのアプリケーションを作成し,uWSGIとnginxで公開しようとした際に,いろいろ詰まったので,どのようにして回避しようとしたのかを共有したいと思います.

前準備

yum install -y nginx
yum install -y https://centos7.iuscommunity.org/ius-release.rpm
yum install -y python36u python36u-libs python36u-devel python36u-pip
yum groupinstall "Development Tools"
yum install -y uwsgi

あとGoogle Domainsなどで適当にドメインを取ってきて,サーバのIPと結びつけておきましょう.

簡素なFlaskアプリを作る

ぼくの環境ではデフォルトでrootとしてLoginしていたので/root/ディレクトリ下にアプリケーション作るとします.

mkdir app
cd app

Python3.6の仮想環境を構築して,適当にアプリケーションを組みます

python3.6 -m venv env
. env/bin/activate
pip3 install flask uwsgi
vim app.py

/root/app/app.py

from flask import Flask

app = Flask(__name__)

@app.route("/")
def index():
    return "Yay! You're on Flask!"

if __name__ == "__main__":
    app.run(host="0.0.0.0" ,port=5000)

とりあえず,ファイアウォールで5000番ポートを開放してhttp://[サーバのIP]:5000 にアクセスします

firewall-cmd  --permanent --add-port=5000/tcp
firewall-cmd --reload
python3 app.py

アクセスして,Yay! You're on Flask!と表示されたら成功です2

Nginx

サーバの設定を行います.ただし,上のようにHTTPプロキシではなくsocketでつなげます.

vim /etc/nginx/conf.d/app.conf

/etc/nginx/conf.d/app.conf

server {
    listen 80;
    server_name [取ってきたドメインの入力];
    charset utf-8;
    client_max_body_size 75M;

    location / { try_files $uri @myapp; }
    location @myapp {
            include uwsgi_params;
            uwsgi_pass unix:///root/app/uwsgi.sock;
    }
}

恥ずかしながら,これを調べる前まではHTTPプロキシでポートを消費して通信する方法しか知りませんでした...

uWSGI

WSGI(Web Server Gateway Interface)*1は、Pythonで書かれたWebアプリケーションとアプリケーションサーバーのあいだのインターフェースを規定した仕様(PEP333)です。WSGI仕様に従ったWebアプリケーションを乗せることができるWebサーバーをWSGIサーバーと呼ぶことにします。

uWSGI(http://projects.unbit.it/uwsgi/)は、WSGIサーバーですが、既存のWebサーバーに機能を追加する形になっています。「WSGIサーバー = Webサーバー + WSGIコンテナ」とするなら、WebサーバーとしてNginx、WSGIコンテナとしてuWSGIを使った組み合わせが Nginx+uWSGI ということになります。

(http://d.hatena.ne.jp/m-hiyama/20120312/1331513519)

とのことです.

(ドキュメンテーションを読む気が起きなくて)この周りの設定で消耗しました...

vim uwsgi.ini

/root/app/uwsgi.ini

[uwsgi]

module = app
socket = /root/app/%n.sock
chmod-socket = 666
callable = app

いろいろ試行錯誤しましたが結果4行で済みました. moduleはapp.pyappの部分,callableはapp.pyのスクリプトの中で呼び出しているapp = Flask(__name__)のappです.スクリプト名等を変更している場合は適宜変更してください.

service nginx restart
adduser centos
passwd centos # パスワードを入力...

rootだと動かないようなのでアプリを動かす用のユーザを作成しました.(もっと良いやり方があるかもしれないけど…)

次に権限とファイアウォール周りを変更します.

chown centos:centos -R .
chmod 555 ..
firewall-cmd  --permanent --add-port=80/tcp
firewall-cmd --reload

起動させてみましょう

uwsgi --ini uwsgi.ini --gid centos --uid centos

URLにアクセスしてエラーが出なければ成功です.

デーモン化させて終わり

最後にこの簡単なプロセスをデーモン化させます.

といっても簡単でuWSGIの起動時に--daemonize /path/to/log/fileと設定するだけでよいようです. ついでに,サービスを止める際にpidで止められるようにpidを/tmp/に出力させます.

[uwsgi]

module = app
socket = /root/app/%n.sock
chmod-socket = 666
callable = app
daemonize = /var/log/uwsgi/%n.log
pidfile = /tmp/myapp.pid
uid = centos
gid = centos

これで準備完了であとは

uwsgi --ini uwsgi,ini

で起動させて

uwsgi --stop /tmp/myapp.pid

で停止させることができます

あ~,俺デーモンになっちゃったよぉ~

日記

研究室の面接が終わり,まああとは天に祈るだけとなり,ある程度緊張も解けてきたので久々に日記っぽい記事を書きます.

自然言語処理

ゼロから作るDeep Learning ? ―自然言語処理編

ゼロから作るDeep Learning ? ―自然言語処理編

自然言語処理に対して,今まで嫌だな~,やりたくないな~と避けてきたのですが,いい加減知識の裾野を言語処理領域まで広げるべきなんじゃないかと,つい最近のBERTだとかちょっと前のAttention/Transformerだとかの興盛の様子を見て購入しました.

1章は前冊の復習であったので,おおよそほとんど既知であったのですが,2章からは

  • 2章  カウントベースの分散表現
  • 3-4章 推論ベースの分散表現
  • 5-6章 RNN,LSTM
  • 7章  RNNによる文章生成
  • 8章  Attention機構

と,昨今のNLP論文を追うために必要な基礎理論を凡そ学習できる良い本だと思う.


  1. ConoHaのVPSを選択した重要な理由として美雲このはちゃんが可愛いという点があります.

  2. お察しの通りYay! You’re on Railsのパクリです

GCPの無料枠で何か作りたかった::10/20

https://i.gyazo.com/5756071422d6a1014e2d0966d6abbc7a.png

こんなものを作りました.

http://instasan.kcs18.net/

追記:開発途中なので容赦なく止まったり不可解な行動を起こす場合があるので優しく扱おう

GCPのプロモが折角あるので利用することにしました.

Github上のmasterにプッシュするとCloud Buildが実行されて,レポジトリ内のDockerfileに基づきコンテナが自動生成されます.そのコンテナをKubernetes Engineで運用しています.ついでにGoogle Domainでドメインを取得したのでGoogle上ですべてが完結してしまっています.恐ろしや.

Kubernetes難しくないっすか?なんもわからんのだが.

日記::10/12

生のJSが割とすごい*1

Web開発は全くの初心者なんで,何を使えば何ができるのか,何かを実現するためには何を使いどのように実装を行えばいいのかといったノウハウが全く分かりません.

インターン(バイト)先でWeb周りの開発補助をすることとなったので,Django等の学習を行っているのですが,MVCアーキテクチャ(のM)が要らない,DBの要らない,割と静的なWebページを経験として作ってみようと思い,FlaskとJavascript,CSS,HTMLという古代Web開発三種の神器を用いてちょっとしたWebアプリ(と呼んでいいのかは知らない)を作成してみました.

止まってます!!!

http://35.231.208.50:40001/

GCP上でFlaskを走らせているだけなので見事なまでの生IPのURLです.こういう感じの公開の仕方がいいのかどうかすらわからん.

Image from Gyazo

MNISTのモデルを使ってcanvasに書いた画像を分類しているだけです.JS+jqueryでこんな感じのことができるのって結構スゴイ気がします.

ページの上に書いてある通り精度はクッソポンコツなので気にしないでね.

Image from Gyazo

MNISTには絶対要りませんが色を選べます.


github.com

コードです.

サーバ側をPythonにさせることで機械学習系Webアプリを簡単に作れるのはうれしいですね.

11月にある三田祭に向けて,なにかサークルの面々と面白そうなアプリを作ってみたいなあという所存です.

*1:こう書くとかなりいかがわしい

Large Scale GAN Training for High Fidelity Natural Image Synthesis(BigGAN)::10/6

https://arxiv.org/pdf/1809.11096.pdf

-1. AbstractのAbstract

緻密に色々試したら,つぎは計算機で殴ろうな!!

0. Abstract

いわゆる,BigGANという名前で,少し話題になっているGANです.DeepMindのインターンの成果物らしく自己アテンションを組み込んだSelf-Attention GANをベースに様々な最適化手法が取り込まれているようです.Appendix Gにあるように「こうしたら悪い結果が出た」というのを掲載しているので実装等の際には非常に有用となるかと思います.チャンネル数とバッチサイズの大小と精度についても述べられており,双方ともに上げた方が良い結果が望めるとのこと(やはり計算資源のデカさこそ正義らしい).潜在変数$z$のサンプリング範囲をあらかじめ決めておいて,その範囲外のサンプルを切り捨てる"Transaction Trick"を使用することで多様性を実現しているようです.Scoreとしては従来のState-of-the-art(SoTA)であるInception Score(IS)とFIDの値52.52,18.65よりも大幅に高精度の166.3,9.6を記録したという驚きの結果が出ています.

1. Introduction

GANの画像生成で往々にして話題になるのは生成画像のそっくり度(Fidelity)と多様性(Variety)です.

以下の3つを考慮することでこれの更新を達成しようと試みたようです.

  1. GANはスケーリングにより大きな影響を受ける.つまるところ,バッチサイズ,パラメータを4倍程度にすることでGANは大幅に利益を受けることを示す.計算資源こそ正義.加えて,スケーラビリティを向上させる簡単なGANアーキテクチャを導入し,また,正規化手法を新たに導入しコンディショニングを改善させることで,大幅なパフォーマンスの向上を図った.
  2. モデルのFidelityとVarietyのトレードオフを細かに制御ができる単純な"Truncation Trick"を導入した.
  3. 大規模なGANを訓練させえる際に特有な不安定性を考察した.経験的にそれらを特徴付けすることで,既存技術と新規技術を融合させることでこの不安定性がある程度軽減できることを立証した.

3. Scaling Up GANs

この節では,GANが大きなバッチサイズとパラメータに対して利益を得るために必要な,GANのトレーニングのための手法を考察するようです.

SA-GAN(Zhang et al.,20181)がベースで,GANの損失関数としてはヒンジ関数(Lim & Ye,20172 )を利用したらしい.$G$に分類の情報を与えるために,class-conditional Batch Norm(Dumoulin et al.,20173)を採用し,$D$にはprojection(宮戸(PFN)&小山(立大),20184 )を採用.Spectral Norm(宮戸,20185)を$G$にも適用(SA-GANと同じ手法)を取ったが,学習率を半分にし,$D$の学習と$G$の学習の回数の比を$2:1$にした.

また,PGGAN(Karras et al.,20186)でも取られている$G$の重みを動的に制御するシステムを減衰率$0.9999$で適用し,重みの初期化は直交初期化(Orthogonal Initialization/Saxe et al.,20147)を使った.

計算資源は昨今話題のTPU v3(Tensor Processing Unit/Google,2018)を使った.さらに512四方の大きな画像でも,PGGANのような逐次的に高解像度へと発展させていく機構が必要ないと言えたらしいです.


まず,Batch Sizeをもとの論文(SA-GAN)の$256$から2,4,8倍へと大きくさせていた結果Batch Sizeに応じてスコアが漸増しており,計算資源こそが正義ということが分かります.Batchが大きいほど入力の多様性が増えるわけですので所謂モード崩壊8が生じにくくなると単純に推測ができるようです.また,Batch Sizeの増加は単純にBatch par Epochが減少するので訓練時間がより短くなるという嬉しい副作用があります.しかしながら,場合によっては学習が失敗してしまうという結果が観測されたようで,上記の結果は学習の崩壊前の最善のウェイトを使った際の結果とのことです.学習の安定化についての考察は第4節にて考察するそうです.

次に各々のレイヤのチャンネル数を150%($64\rightarrow 96$)にしてパラメータ数を大体倍にした際もスコアの上昇がみられたようです.

また,$G$で使われている条件付きBatch Normは多くのパラメータを含んでいる.そこで,各レイヤで埋め込みを行う(宮戸,2018[^4])のではなく,各々のレイヤのバイアスとゲインを線形投影させる,共有埋め込み(Shared Embedding)を用いてパラメータと学習時間(=計算コスト)を削減させた.

さらに,はじめの層以外にノイズを入力させる技法を採用した.この目的としては潜在空間を利用して異なる解像度や異なる階層のレベルで影響を及ぼすことが目標らしい9.

3.1. "Truncation Trick"によるVarietyとFidelityのトレードオフ

GANは入力潜在変数分布$p(z)$を任意の分布に設定することができるが,ほとんどのGANの論文ではたいてい正規分布$\mathcal{N}(0,1)$や,一様分布$\mathcal{U}[-1,1]$などからノイズを生成している.この事前分布の選定によるパフォーマンスの違いをオマケEにてまとめてあるようだ.

数々の事前分布で実験を行ったところ,驚くべきことに,学習時とサンプリング(予測)時で異なる分布からサンプリングしたほうが良い結果となることが分かったようだ.学習時には$z\sim\mathcal{N}(0,1)$で,生成時には打切り正規分布(Truncated Normal:最頻値からの特定範囲外のサンプルを棄却した正規分布)を用いた際に最も良い結果となることが観測された.このテクニックをTruncation Trickと呼ぶことにする.

打切りの閾値は生成画像の多様性とそれっぽさに大きく影響を与える.上の図は閾値を$2,1.5,1,0.5,0.04$10とした場合である.

打切りの閾値を操作することで,事後的かつ詳細に(Fine-Grained Post-Hoc)生成画像の多様性とそれっぽさを変えることが可能となる.

実際にこのトレードオフを可視化すべくISとFIDという評価指標を用いて実験を行ったようだ.

ISはInceptionモデルで識別しやすい画像でかつ,識別されるラベルの多様性が豊富なほど11高スコアを記録するものであり,FIDはGANの生成する画像の品質を評価する12もので,真の分布の画像集合と再現分布の画像集合間の距離を評価するものであるので,低スコアである方が望ましい.

(実際に評価した図はオマケDの図16に載っている)

しかし,大きなモデルの中には打切りがうまく作用せずサチった(Saturated)画像を生成してしまう場合があるようだ(Saturation Artifact).これに対応するために,$G$を平滑化させるように条件付けし,切り捨てに対する順応性を上げる必要性,つまりは,潜在変数$z$の全空間に置いて良好な出力に射影されるようにさせる.これを達成すべく,直交正則化(Orthogonal Regularization,Brock et al,2017[^6])を適用した: $$ R_{\beta}(W)=\beta||W^{T}W-I||^{2}_{F} $$ ただし,$||\cdot||_{F}$は行列の要素ごとの平方和の平方根であるフロベニウスノルムで,$W$は重みパラメータ,$\beta$はハイパーパラメータである.この正則化は制限が強すぎる[^4]ことが指摘されており,そのため,より緩い制約となるような正則化項を探索したところ次のような正則化項が最も有用であることがわかった. $$ R_{\beta}(W)=\beta||W^{T}W\odot({\bf 1}-I)||_{F}^{2} $$ ただし,$\odot$は行列のアダマール積(要素ごとの積)で,${\bf 1}$は全要素が1である行列である.

$\beta=10^{-4}$程度の軽い正則化を課したところ,パフォーマンスの向上が見受けられた.

4. 分析

4.1. Characterizing Instability: Generator

GANの分析の際は簡単な分布を用いて考察を行うことが多いが,ここでは大規模の問題を直接解析したいので,Odena et al(2018)13 でとられているように,学習中の重みと勾配,損失を監視して学習の崩壊する要因となりうる予兆となる指標を探すことにする.

観測の結果各重み行列の第三位までの特異値$\sigma_0,\sigma_1,\sigma_2$が,崩壊の観測に対して有用であることが明らかとなった.行列の特異値は宮戸 et al.(2018)[^4]で使われている冪乗法を拡張したAlrnoldiの反復法(Golub et al.,200014 )を用いて計算できる.

実際の学習の際にこの(Spectral Normを適用する前の)特異値の変化をプロットしたところ学習崩壊の際に特異値が急に爆発することが観測された.

この病的なふるまいが,崩壊の原因なのか症状なのかを究明するため,$G$に対してスペクトル爆発を明示的に抑制させる為の追加の条件を課すことにしたようだ.

スペクトル爆発を制御する方法としては,重み行列$W$をとして,第一特異ベクトルと特異値$u_0,v_0,\sigma_0$として,$\sigma_0$に対してクランプされる値を$\sigma_{clamp}$とすると,重みを次のようにして更新する. $$ W=W-\max(0,\sigma_0-\sigma_{clamp})v_0u_{0}^{T} $$ ここでクランプ値は次のようにして設定することにする.

  1. 単に固定値$\sigma_{reg}$を設定してクランプする
  2. 第二特異値$\sigma_1$に対する比$r\cdot sg(\sigma_1)$.ただし,$sg(\cdot)$は正則化により$\sigma_1$を増加させないようにするための勾配に対する停止操作である15

クランピングである程度のパフォーマンス改善や,爆発防止が見込める場合があるが,完全に学習の崩壊の防止を行うことはできないようだ.

スペクトルノルムは学習安定化のための重要な要素であるが,現在のテクニックではうまく安定化させることは難しいようだ(悲しい).

4.1の補足: 特異値

特異値分解は,正方行列にのみ適応可能な固有値分解を縦横が一致しない一般の矩形行列にたいして適用できるように拡張させたものであるといえます.

特異値分解は行列$M\in K^{m\times n}$(体$K$は実数や複素数を考える)をユニタリ行列$U\in K^{m\times m},V\in K^{n\times n}$と,半正定値行列$MM^{\star}$(ここで$A^{\star}$は$A$の転置+複素共役行列(随伴行列ともいう))の正の固有値の平方根$\sigma_1\geq\sigma_2\geq...\geq\sigma_r>0$が存在して,$q=\min(m,n),\sigma_{r+1}=...=\sigma_{q}=0$と置いて,$\Delta=\text{diag}(\sigma_1,\sigma_2,...,\sigma_q)$として,$\Sigma\in K^{m\times n}$を左上に$\Delta$を置いて,残りの部分を0埋めした行列として, $$ M=U\Sigma V^{\star} $$ という形に分解することをいい,$\sigma_r$は特異値と呼ばれます.

固有ベクトルが正規直交基底となったことと同様に特異値分解でも行列$V$が入力空間における正規直交基底となり,$U$が出力空間における正規直交基底と等しいといえ,$\Sigma$は成分のゲインであると考えることができる.基本的に特異値は降順に並べることが多く,この場合$\Sigma$が一意に定まることが知られています.

上記特性から特異値は $$ \begin{eqnarray} Mv&=&\sigma u\ M^{\star}u&=&\sigma v \end{eqnarray} $$ が成り立つ$u,v$組に対応する$\sigma$であるといえ,$u$を左特異ベクトル,$v$を右特異ベクトルと呼ばれます.

特異値分解を行うメソッドはnumpyに実装されておりnp.linalg.svg(Matrix)を実行することで返り値として$U,\Sigma,V^{\star}$を得ることができます.

たとえば, $$ A=\left[ \begin{array}{ccc} -\sqrt{2}^{-1} & 1 & -1\\ -\sqrt{2}^{-1} & -1 & 1 \end{array} \right] $$ の特異値は $$ \Sigma=\text{diag}(2,2) $$ のように表すことができ,

import numpy as np
A = np.array([[-(2)**(.5), 1, -1], [-(2)**(.5), -1, 1]])
U, s, V = np.linalg.svd(A)

とすることで,実際に特異値分解を行うことができます.

print(U)
print(s)
print(V)
# [[-1.00000000e+00 -2.22044605e-16]
#  [-2.22044605e-16  1.00000000e+00]]
# [2. 2.]
# [[ 7.07106781e-01 -5.00000000e-01  5.00000000e-01]
#  [-7.07106781e-01 -5.00000000e-01  5.00000000e-01]
#  [-1.75725255e-16  7.07106781e-01  7.07106781e-01]]

さて,特異値を実際に近似で求めるアルゴリズムの中で,Spectral Normの論文で紹介されているものとして冪乗法(Power Method)があります.冪乗法は以下からなるアルゴリズムです


アルゴリズム 1:冪乗法


入力:$A\in K^{m\times n}$と反復回数$N$

出力:最大(第一)特異値$s$16 ,とその左及び右の特異ベクトル$u,v$

let $q\in K^{n}, \text{where,}~q_i \leftarrow \mathcal{N}(0,1),||q||_2=1$

for $i$ in $0..N$:

​$~~~~~~~~q\leftarrow A^{T}Aq$

$v \leftarrow q/||q||_2$

$s\leftarrow ||Av||_2$

$u\leftarrow Av/s$

return $u,s,v$


これをPythonで実装させると,

def power_method(A, num_iter=1):
    q = np.random.normal(size=A.shape[1])
    q = q/np.linalg.norm(q)
    for i in range(num_iter):
        q = A.T.dot(A).dot(q)
    v = q / np.linalg.norm(q)
    s = np.linalg.norm(A.dot(v))
    u = A.dot(v)/s
    return u, s, v

このようになり17,実際に先ほどの例を実行させると,

print(power_method(A, 5))
# output:
# (array([ 0.6678161 , -0.74432631]),
# 1.9912766893959086,
# array([0.07617649, 0.70299157, 0.07617649, 0.70299157]))

と,およそnumpyの実行結果と一致していることが確認できます.

対して,この論文で紹介されているArnoldi法の実装はコレ18が詳しいので同様に実装することができます.

4.2. Characterizing Instability: The Discriminator

$G$と同様に$D$に対してもスペクトルノルムの分析を行い,制約を課して学習を安定させることを考える.

$D$の特異値は$G$に比較して,ノイズが多い傾向にあり,また崩壊時は$G$と違って値が爆発的に増加するのではなく漸増量が増加する傾向にある.

ところどころ$D$の特異値に観測される尖頭値についてまず考察をする.この尖頭値から,周期的に非常に大きな勾配を受け取るようなことを示唆しているように思えるが,$D$のフロベニウスノルムが滑らかであることを考えると,この尖頭値現象はいくつかの主位特異ベクトル方向に集中しているといえる.つまり,尖頭値では$G$が$D$を激しく攪乱させるような画像を生成されたが故の敵対的プロセス下での最適化の結果によるものだと仮定できる.

このスペクトルノイズが不安定性と関連があるかどうか調べるために,$D$のヤコビアンを明示的に制限するような勾配罰則を適用して正則化する.この実験ではMescheder et al.(2018)19の,$R_1$ゼロ中心勾配罰則(Zero-centered Gradient Penalty)を導入した.

$$ R_1:=\frac{\gamma}{2}\mathbb{E}_{p_{\mathcal{D}}(x)}[||\nabla D(x)||_F^{2}] $$

$\gamma$を$10$とすることで学習が安定化し,スペクトルの平滑さが向上したがパフォーマンスは著しく低下した.特にIS(多様性)の方のスコアが45%ほど落ち込む結果となった.むろん$\gamma$を緩和することでパフォーマンスは向上するが次第にスペクトルの尖頭が目立つようになり始め,$\gamma=1$でもISは20%のパフォーマンス低下が生じるという結果となった.

同様の実験を直交正規化[^6],DropOut20,L2正則化で行ったが同様の結果となった.

5. 実験

ImageNetとJET-300Mで実験を行ったようです.

ImageNet

結果が左から

  • Truncation Trickを行わなかった場合のスコア
  • FIDを最小化したときのスコア(最大精度)
  • ISの検証セットを利用したときのスコア(汎化の客観性)
  • ISを最大化したときのスコア(最大多様性)

となっている.

JET-300M

256四方の場合における実験結果のようです.


オマケにはパラメータ等調整の際の試行錯誤の結果が載っており再現実験等をする場合非常にありがたいかもしれませんが,512四方の画像をチャネル96,バッチ数2048で実行できる計算機でぶん殴る必要があるでしょう.つらい.


  1. Self-Attention Generative Adversarial Networks https://arxiv.org/abs/1805.08318

  2. Geometric GAN,https://arxiv.org/abs/1705.02894

  3. A Learned Representation For Artistic Style,https://arxiv.org/pdf/1610.07629.pdf

  4. cGANs with Projection Discriminator,https://arxiv.org/abs/1802.05637

  5. Spectral Normalization for Generative Adversarial Networks,https://arxiv.org/abs/1802.05957

  6. Progressive Growing of GANs for Improved Quality, Stability, and Variation,https://arxiv.org/abs/1710.10196

  7. Exact solutions to the nonlinear dynamics of learning in deep linear neural networks,https://arxiv.org/abs/1312.6120

  8. $G$が最頻(mode)単峰分布のみをめちゃくちゃよい精度で生成するように学習して多様性のある画像が生成されなくなる現象

  9. ここがいまいち理解できなかった.初期層以外にノイズを入力することによる利点が何かを別途別文献参照して考察したい

  10. 図の個数と閾値の数が合ってない?

  11. http://bluewidz.blogspot.com/2017/12/inception-score.htmlこの記事が詳しい

  12. http://bluewidz.blogspot.com/2017/12/frechet-inception-distance.html同じくこの記事が詳しい

  13. Is Generator Conditioning Causally Related to GAN Performance?https://arxiv.org/abs/1802.08768

  14. Journal of Computational and Applied Mathematics,(検索結果がArnoldiしか出ないが原文ママである,http://shnya.hatenablog.jp/entry/2016/03/12/100507が詳しい)

  15. どういう操作をすればいいのかいまいちわかりませんでした…

  16. 行列の最大の特異値のことは行列のスペクトルノルム(Spectral Norm)といわれるものです.行列に対する$L_2$ノルムであり,$||A||_{2}=(Maximum~Singular~Value)=\max_{||x||_{2}\neq0}||Ax||_2/||x||_2$であることが示されています.特異値が入力に対するゲインと考えるとSpectral Normは入力を極端に大きくスケーリングしないように調整している.とでもいえるのでしょうか?(←この私の考え方が正しいかどうかは判断しかねますが)

  17. https://blog.tsurubee.tech/entry/2017/11/05/225803を参照にした.

  18. http://ry0u.github.io/post/2018-06-16-arnoldi-and-eigenvalues/

  19. Which Training Methods for GANs do actually Converge? https://arxiv.org/abs/1801.04406

  20. Dropout: A Simple Way to Prevent Neural Networks from Overfitting http://jmlr.org/papers/volume15/srivastava14a.old/srivastava14a.pdf

日記(9月)

日記のレイアウトを大幅に変えました.しっくり来たのでこれで行こうかと思います.

とあるところでバイト(インターン)を始めました.実務に触れる機会としてはインターンは本当に有用ですよね.実際に何か学問を学ぶという場合でも,先に自分の求める分野のエキスパートがどんな道具を使ってどんな成果物を出しているかとかをしっかり確認することの重要性を学んだ気がします.

Web開発

Image from Gyazo

実用性とかパフォーマンスはともかく簡単にアニメーションを作ることができてホント脱帽という感じでした.たのしい.

ウェブ開発は闇が深そうだと勝手に思ってました.しかし,ある程度触ってみてやりたいことを調べたら98%ライブラリがあって先人たちの知恵って偉大だなあと思うと同時に闇が少し晴れた気がします.相変わらず靄程度の闇はありますけど.

  • ReactJS
  • Ruby on Rails(Progate)
  • Django(業務で使うため)

言語もばらばらですけど所謂RESTだとかMVCだとかの思想はどれも同様であるので,割と第二言語以降の習熟は外国語と違って*1早くできるなあと思っています.文法記法とAPIに依存した思念さえ習得してしまえば案外簡単にできてしまうものなんですね.

Python Django 超入門

Python Django 超入門

読んでいます.

prog-8.com

やっています.

競技プログラミング

基本的なアルゴリズムだとか,柔軟な思考力とかを身に着けたくてAtCoderと蟻本を読み始めました.

プログラミングコンテストチャレンジブック [第2版] ?問題解決のアルゴリズム活用力とコーディングテクニックを鍛える?

プログラミングコンテストチャレンジブック [第2版] ?問題解決のアルゴリズム活用力とコーディングテクニックを鍛える?

beta.atcoder.jp

Rating: 156 -> 414

せめて水色コーダーくらいまでには今年中に到達したいなあと思っています.(青はさすがに厳しい険しい)

練習がてら所謂グラフ理論の「グラフ」を扱うクラスを作ってみたのですが木の半径を求めるのがどうも遅くて$N=100000$で8秒くらいかかってしまう要改善なコードですがDドライブの中で供養するのはかわいそうだったのでGistに眠らせることにしました.

gist.github.com

慶応義塾大学のココがすごい

9/24より後期の講義が始まったらしく,矢上の坂登りの日々が始まりました.

授業のひとつに簡易的なロボコンのようなものがあって,今までハード分野に非常に疎かったので触るいい機会となりました*2

他には最適化数理という名のアルゴリズムの講義やシステム計測という名の待ち行列理論講義なんかがあり結構面白いです.

後期は実験がなくて最高だぜ!

数学

論文を読むにあたり基本的な数理が全く分かってないなと感じ,三年の秋にしてようやく数学に手を出しました.

集合と位相 (数学シリーズ)

集合と位相 (数学シリーズ)

さしあたって集合と位相を読んでいます.「Zornの補題」「選択公理」「整列可能定理」が同値ってスゴイっすね!何言ってっか全く理解できねぇけど!

機械学習

pix2pixいまだにうまくいかない...パラメータチューニングでおわるお前の人生...

ブログ

3回しか書いてない.日記とは何だったのか.

fgjiutx.hatenablog.com

fgjiutx.hatenablog.com

fgjiutx.hatenablog.com

*1:関連するどうでもいい余談として講義でロシア語を取ったのですができる気がしない.3とか|-Oとか文字の時点で挫折しそう

*2:Q.ロ技研は…? A.聞かないで

初心に帰って変分ベイズ(2)::9/17

fgjiutx.hatenablog.com この記事の続きです

4.2. 変分事後分布の算出

モデルを変分ベイズ法で扱うために,まず,全確率変数の同時分布を求めます. $$ p({\bf X},{\bf Z},{\bf \boldsymbol\pi},{\bf\boldsymbol \mu},{\bf \Lambda})= p({\bf X}|{\bf Z},{\bf \boldsymbol\mu},{\bf \Lambda})p({\bf Z}|{\bf \boldsymbol\pi})p({\bf\boldsymbol \pi}) p({\bf\boldsymbol \mu}|{\bf \Lambda})p({\bf \Lambda}) $$ 観測変数は${\bf X}$のみであるので,変分近似は,潜在変数とパラメータに分解して,以下のように考えられます. $$ q({\bf Z},{\bf\boldsymbol \pi},{\bf \boldsymbol\mu},{\bf \Lambda})= q({\bf Z})q({\bf \boldsymbol\pi},{\bf \boldsymbol\mu},{\bf \Lambda}) $$ 前述のように,因子$q({\bf Z})$は,他方の因子の期待値を求めることで,最適化でき, $$ \ln q^{\star}({\bf Z})=\mathbb E_{\pi,\mu,{\bf \Lambda}} [\ln p({\bf X},{\bf Z},{\bf \boldsymbol\pi},\boldsymbol\mu,{\bf \Lambda})]+\text{const} $$ 同次分布の分解式から,${\bf Z}$にかかわる項だけを抜き出すことで,以下のように式変形をすることができます. $$ \ln q^{\star}({\bf Z})=\mathbb E_{\boldsymbol\pi,\boldsymbol\mu,{\bf \Lambda}} [\ln p({\bf Z}|\boldsymbol\pi)]+ \mathbb E_{\boldsymbol\pi,\boldsymbol\mu,{\bf \Lambda}} [\ln p({\bf X}|{\bf Z},\boldsymbol\mu,{\bf \Lambda})]+\text{const} $$ ここに,実際の式を代入することで, $$ \ln q^{\star}({\bf Z})=\sum_{n=1}^{N}\sum_{k=1}^{K}z_{nk}\ln \rho_{nk} + \text{const} $$ となります.ただし, $$ \ln\rho_{nk}=\mathbb E[\ln \pi_k]+\frac{1}{2}\mathbb E[\ln {\bf |\Lambda_k|}] -\frac{D}{2}\ln 2\pi-\frac{1}{2}\mathbb E_{\boldsymbol\mu_k,{\bf \Lambda}_k} [({\bf x}_n-\boldsymbol\mu_k)^{T}{\bf \Lambda}_k({\bf x}_n-\boldsymbol\mu_k)] $$ で,$D$はデータの次元です.

式の指数を取ると, $$ q^{\star}({\bf Z})\propto\prod_{n=1}^{N}\prod_{k=1}^{K}\rho_{nk}^{z_{nk}} $$ となり,この分布は正規化されている必要があるので, $$ r_{nk}:=\frac{\rho_{nk}}{\displaystyle \sum_{j=1}^{K}\rho_{nj}} $$ とすることで, $$ q^{\star}({\bf Z})=\prod_{n=1}^{N}\prod_{k=1}^{K}r_{nk}^{z_{nk}} $$ が,得られます.この最適化の式は,式中に他のパラメータが出てきているため,相互依存関係にあり双方を繰り返し更新していく必要があります.

次に,$q^{\star}({\bf \pi},{\bf \mu},{\bf \Lambda})$を考えます.

$q^{\star}({\bf Z})$と同様,一般的な変分推定の式$\ln q^{\star}(\boldsymbol\pi,\boldsymbol\mu,{\bf \Lambda})=\mathbb E_{{\bf Z}}[\ln p({\bf X},{\bf Z},{\bf \boldsymbol\pi},\boldsymbol\mu,{\bf \Lambda})]+\text{const}$から, $$ \ln q({\bf \boldsymbol\pi},{\bf \boldsymbol\mu},{\bf \Lambda})= \ln p(\boldsymbol\pi)+\sum_{k=1}^{K}\ln p(\boldsymbol\mu_k,{\bf \Lambda}_k)+\mathbb E_{{\bf Z}} [\ln p({\bf Z}|\boldsymbol\pi)]+\sum_{k=1}^{K}\sum_{n=1}^{N}\mathbb E[z_{nk}] \ln \mathcal N({\bf x}_n|\mu_k,{\bf \Lambda}_k^{-1})+\text{const} $$ と分解されます.

分解された式は,$\pi$を含む項と,$\boldsymbol\mu,{\bf \Lambda}$を含む項に分割することができるので,$q({\bf \boldsymbol\pi},{\bf \boldsymbol\mu},{\bf \Lambda})$は, $$ q({\bf \boldsymbol\pi},{\bf \boldsymbol\mu},{\bf \Lambda})=q(\boldsymbol\pi)\prod_{k=1}^{K}q(\boldsymbol\mu_k,{\bf \Lambda}_k) $$ となり,$q^{\star}(\pi)$は $$ \ln q^{\star}(\boldsymbol\pi)=\ln p(\boldsymbol\pi)+\mathbb E_{\bf Z}[\ln p({\bf Z}|\boldsymbol\pi)]+\text{const} $$ から,$p({\bf \pi})$は, $$ \begin{eqnarray} \ln p({\bf \boldsymbol\pi})&=& \ln \left(C(\alpha_0)\prod_{k=1}^{K}\pi_{k}^{\alpha_{0}-1}\right)\\ &=&(\alpha_0-1)\sum_{k=1}^{K}\ln \pi_k+\text{const} \end{eqnarray} $$ となって,$\mathbb E_{\bf Z}[\ln p({\bf Z}|\boldsymbol\pi)]$は $$ \begin{eqnarray} \mathbb E_{\bf Z}[\ln p({\bf Z}|\boldsymbol\pi)]&=&\mathbb E_{{\bf Z}} \left[\ln\left(\prod_{n=1}^{N}\prod_{k=1}^{N}\pi_{k}^{z_{nk}}\right)\right]\\ &=&\sum_{k=1}^{N}\sum_{k=1}^{N}\mathbb E_{{\bf Z}}[z_{nk}]\ln \pi_k+\text{const} \end{eqnarray} $$ $z_{nk}$は0もしくは1をとるので期待値は容易に計算でき,$\mathbb E_{\bf Z}[z_{nk}]=r_{nk}$となるので, $$ \mathbb E_{\bf Z}[\ln p({\bf Z}|\boldsymbol\pi)]=\sum_{k=1}^{N}\sum_{k=1}^{N}r_{nk}\ln \pi_k+\text{const} $$ となるので,結局, $$ \ln q^{\star}(\boldsymbol\pi)=(\alpha_0-1)\sum_{k=1}^{K}\ln \pi_k+\sum_{k=1}^{N}\sum_{k=1}^{N}r_{nk}\ln \pi_k+\text{const} $$ となります.ここで,指数を取ると, $$ q^{\star}(\pi)=\prod_{k=1}^{K}\pi_{k}^{\alpha_0+\sum_{n=1}^{n}r_{nk}-1} $$ より,$N_k=\sum_{n=1}^{n}r_{nk}$,${\bf \alpha}_k=\alpha_0+N_k$とすると,$q^{\star}(\pi)$はディリクレ分布に従うことが分かります. $$ \ln q^{\star}(\boldsymbol\pi)=\text{Dir}(\boldsymbol\pi|\boldsymbol \alpha) $$ つぎに,$q(\boldsymbol\mu_k,{\bf \Lambda}_k)$は,乗法定理から$q^{\star}(\boldsymbol \mu_k,\boldsymbol \Lambda_k)=q^{\star}(\boldsymbol\mu_k|\boldsymbol\Lambda_k)q^{\star}(\boldsymbol\Lambda_k)$となるので,変分推定式から$\boldsymbol \mu_k,\boldsymbol \Lambda_k$の項を抜き出して(長い)式変形をすると,ガウス-ウィシャート分布となることが分かります.

ここで,ガウス-ウィシャート分布を考えるうえで次のコンポーネント(クラスタ)ごとの加重平均した統計量を用いると式の簡便化が見込めます. $$ \begin{eqnarray} \bar{{\bf x}_k}&=&\frac{1}{N_k}\sum_{n=1}^{N}r_{nk}{\bf x}_n\\ {\bf S}_k&=&\frac{1}{N_k}\sum_{n=1}^{N}r_{nk}({\bf x}_n-\bar{\bf x}_k) ({\bf x}_n-\bar{\bf x}_k)^{T} \end{eqnarray} $$ この補助的な統計量を用いると, $$ q^{\star}(\boldsymbol\mu_k,\boldsymbol\Lambda_k)=\mathcal N(\boldsymbol\mu_k|{\bf m}_k,(\beta_k\boldsymbol\Lambda_k)^{-1}) \mathcal W(\boldsymbol\Lambda_k|{\bf W}_k,\nu_k) $$ と表されます.ただし, $$ \begin{eqnarray} \beta_k&=&\beta_0+N_k\\ {\bf m}_k&=&\frac{1}{\beta_k}(\beta_0{\bf m}_0+N_k\bar{\bf x}_k)\\ {\bf W}^{-1}&=&{\bf W}_0^{-1}+N_k{\bf S}_k+\frac{\beta_0N_k}{\beta_0+N_k} (\bar{\bf x}_k-{\bf m}_0)(\bar{\bf x}_k-{\bf m}_0)^{T}\\ \nu_k&=&\nu_0+N_k \end{eqnarray} $$ とします.

これらを振り返ると,

  1. モデルの現在の事後分布から負荷率に当たる$r_{nk}$を計算する
  2. その計算された負荷率をもとに事後分布を更新する

という,最尤推定におけるEMアルゴリズムに似た方法で,最適化を行えばよいということがわかります.また,$r(\propto \rho)$を計算するにあたり,パラメータの分布の期待値が必要となるがこれらは容易に計算することができて, $$ \begin{eqnarray} \mathbb E_{\boldsymbol\mu_k,{\bf \Lambda}_k} [({\bf x}_n-\boldsymbol\mu_k)^{T}{\bf \Lambda}_k({\bf x}_n-\boldsymbol\mu_k)] &=&D\beta_k^{-1}+\nu_k(\bar{\bf x}_k-{\bf m}_0)^{T}{\bf W}_k (\bar{\bf x}_k-{\bf m}_0)\\ \ln\tilde\Lambda_k:=\mathbb E[\ln|\boldsymbol\Lambda_k|]&=&\sum_{i=1}^{D} \psi\left(\frac{\nu_k+1-i}{2}\right)+D\ln2+\ln|{\bf W}_k|\\ \ln\tilde\pi_k:=\mathbb E[\ln \pi_k]&=&\psi(\alpha_k)-\psi(\hat\alpha),~ \text{where}~\hat\alpha=\sum_{k=1}^{K}\alpha_k \end{eqnarray} $$ となります.ただし,$\psi(\cdot)$はディガンマ関数と呼ばれるもので, $$ \psi(a):=\frac{d}{da}\ln\Gamma(a) $$ と定義されます.

これらの結果から, $$ r_{nk}\propto \tilde\pi_k\tilde\Lambda^{1/2}\exp\left(-\frac{D}{2\beta_k} -\frac{\nu_k}{2}({\bf x}_n-{\bf m}_k){\bf W}_k({\bf x}_n-{\bf m}_k)^{T}\right) $$ となります.

これにより,更新を繰り返してモデルを最適化する準備ができました.

4.3. 予測分布の計算

モデルに対してパラメータを最適化させた後は,未知の観測値に対して,予測分布を算出することが目的である場合がおおいので,新たな観測値$\hat x$に対して,分布$p(\hat x|{\bf X})$を求めることを考えます.

上図のようなグラフィカルモデルより考えると,分布の式が分かりやすくなり, $$ p(\hat x|{\bf X})=\sum_{\hat{\bf z}}\iiint p(\hat{\bf x}|\hat{\bf z},\boldsymbol\mu,\boldsymbol\Lambda) p(\hat{\bf z}|\boldsymbol\pi)p(\boldsymbol\pi,\boldsymbol\mu,\boldsymbol\Lambda|{\bf X})\text{d}\boldsymbol\pi \text{d}\boldsymbol\mu\text{d}\boldsymbol\Lambda $$ のように周辺化を行うことで達成されます.

いま,変分近似を用いて実パラメータ分布を$p(\boldsymbol\pi,\boldsymbol\mu,\boldsymbol\Lambda|{\bf X})\approx q(\boldsymbol\pi)q(\boldsymbol\mu,\boldsymbol\Lambda)$と近似したので, $$ \begin{eqnarray} p(\hat{\bf z}|{\bf \boldsymbol\pi})&=&\prod_{k=1}^{K}\pi_k^{\hat z_k}\\ p(\hat{\bf x}|\hat{\bf z},{\bf \boldsymbol\mu},{\bf \Lambda})&=&\prod_{k=1}^{K}\mathcal N({\bf x}_n|{\bf \boldsymbol\mu}_k,{\bf \Lambda}_k^{-1})^{\hat z_{k}}\\ \sum_{\hat{\bf z}} p(\hat{\bf z}|{\bf \boldsymbol\pi})p(\hat{\bf x}|\hat{\bf z},{\bf \boldsymbol\mu},{\bf \Lambda})&=&\sum_{k=1}^{K}\pi_k \mathcal N({\bf x}_n|{\bf \boldsymbol\mu}_k,{\bf \Lambda}_k^{-1})\\ p(\hat x|{\bf X})&\approx&\sum_{k=1}^{K}\iiint\pi_k \mathcal N({\bf x}_n|{\bf \boldsymbol\mu}_k,{\bf \Lambda}_k^{-1})q(\boldsymbol\pi)q(\boldsymbol\mu_k,\boldsymbol\Lambda_k)\text{d}\boldsymbol\pi \text{d}\boldsymbol\mu_k\text{d}\boldsymbol\Lambda_k \end{eqnarray} $$ となり,この積分は,なんとびっくり解析的に解くことができて,混合[多変量スチューデントt-分布](https://en.wikipedia.org/wiki/Multivariate_t-distribution)となります. $$ p(\hat x|{\bf X})\approx\frac{1}{\hat\alpha}\sum_{k=1}^{K} \text{St}(\hat{\bf x}|{\bf m}_k,{\bf L}_k,\nu_k-1+D) $$ ただし, $$ {\bf L}_k=\frac{(\nu_k+1-D)\beta_k}{1+\beta_k}{\bf W}_k $$ です.

4.4. 実装

import numpy as np
from numpy.linalg import slogdet
from sklearn import datasets
from scipy.stats import multivariate_normal, wishart
from scipy.special import gamma, digamma, loggamma
import scipy.stats as stats
import matplotlib.cm as cm
import matplotlib.pyplot as plt
from IPython.display import clear_output
import os

高度な確率分布が多く出てくるため,scipyを多用します.ただし,多変量スチューデントt-分布はscipyにはないので,こちらを参照に自力で実装します.

def multivariate_student_t(x, mu, Sigma, df):    
    x = np.atleast_2d(x) 
    nD = Sigma.shape[0] 

    numerator = gamma(1.0 * (nD + df) / 2.0)
    denominator = (
            gamma(1.0 * df / 2.0) * 
            np.power(df * np.pi, 1.0 * nD / 2.0) *  
            np.power(np.linalg.det(Sigma), 1.0 / 2.0) * 
            np.power(
                1.0 + (1.0 / df) *
                np.diagonal(
                    np.dot( np.dot(x - mu, np.linalg.inv(Sigma)), (x - mu).T)
                ), 
                1.0 * (nD + df) / 2.0
                )
            )
    return 1.0 * numerator / denominator 

さて,ここでこれから3次元テンソルの積を含むような計算を行っていくのですが,numpyで3次元以上の配列(テンソル)の積を扱いたいときに,非常に便利なものとして,numpy.einsumがあります.

その前に...

4.4.1. アインシュタインの縮約記法

たとえば,ベクトル$\boldsymbol a=(a_1,a_2),\boldsymbol b=(b_1,b_2)$の内積は, $$ \boldsymbol a\cdot\boldsymbol b=a_1b_1+a_2b_2 $$ で,一般の$n$次元ベクトルの場合には$\sum$記号を使って, $$ \boldsymbol a\cdot\boldsymbol b=\sum_{i=1}^{n}a_ib_i $$ のように表すことができます.

アインシュタインの記法は大胆にもこの$\sum$を省略して, $$ \boldsymbol a\cdot\boldsymbol b=a_ib_i $$ だとか$a_\mu b^{\mu}$のように書き表す方法です.

np.einsumはこの縮約記法を用いてテンソル積を定義するものです.

例えば,この内積は,

a = np.arange(3) # [0 1 2]
b = np.arange(3) # [0 1 2]

print(np.einsum("i, i ->", a, b)) # 5

と書くことができます.

第一引数がアインシュタインの縮約記法にあたり,->の左側にカンマ区切りでテンソルの縮約記法を並べ,右側に出力テンソルの記法を書きます.ただし,スカラー値の場合は省略します.

これを用いれば,様々な種類のテンソルの乗算を分かりやすく書き下すことができます.

スカラー値×ベクトル
a = 4
b = np.arange(3)

print(np.einsum(", i ->i", a, b)) # [0 4 8]
2次元行列×ベクトル=ベクトル($Ax=b$)
a = np.arange(6).reshape((2, 3)) # [[0 1 2],[3 4 5]]
b = np.arange(3)

print(np.einsum("ij, j -> i", a, b)) # [ 5 14]
$\boldsymbol a^{T}\boldsymbol W\boldsymbol b,\boldsymbol a\in\mathbb R^{k\times j},\boldsymbol b\in\mathbb R^{k\times i},\boldsymbol W\in\mathbb R^{k\times i\times j}$

a = np.arange(6).reshape((2, 3))
W = np.arange(24).reshape((2, 4, 3))
b = np.arange(8).reshape((2, 4))

print(np.einsum("kj, kij, ki -> k", a, W, b)) # [ 156 4844]

4.4.2. クラスの実装

class VariationalMixtureGaussian(object):
    def __init__(self, num_clusters=3):
        '''
        クラスの初期化
        することはせいぜいクラスタ数の指定ぐらい
        '''
        
    def expectation(self, X):
        '''
        X: データ
        期待値を求める過程
        返り値は負荷率 r
        '''
        
    def is_converged(self):
        '''
        収束の判定
        返り値 bool値
        '''
        
    def maximization(self, X, r):
        '''
        X: データ
        r: 負荷率
        事後分布の最大化過程
        '''
        
    def elbo(self):
        '''
        (未解説)変分下限を求める
        返り値 スカラー値
        '''
        
    def fit_init(self, X):
        '''
        X: データ
        最適化の初期化を行う
        '''
        
    def fit_1(self, X, i):
        '''
        X: データ
        i: 最適化を回した回数(printの為に使用)
        1回のみ最適化操作を行う
        '''
        
    def fit(self, X, num_iter=1000):
        '''
        X: データ
        num_iter: 最適化回数の上限
        最適化操作を行う
        '''
        
    def predict(self, X):
        '''
        X: データ
        新たな入力データに対する予測分布の確率を返す
        '''
        
    def classify(self, X):
        '''
        データの分類を行う
        方法は
        1. E過程をしてデータ点の負荷率を計算
        2. argmaxを取って分類する
        '''
        
def multivariate_student_t(x, mu, Sigma, df):
    def __init__(self,num_clusters=3):
        self._param = None
        self.num_clusters = num_clusters

    def expectation(self, X):
        pi = digamma(self.alpha) - digamma(self.alpha.sum())
        det_Wk = np.linalg.slogdet(self.W)[1]
        Lambda = digamma((self.nu - np.arange(2)[:,np.newaxis])/2).sum(axis=0)+\
                 self.ndim*np.log(2) + det_Wk
            
        pi = np.exp(pi)
        Lambda = np.sqrt(np.exp(Lambda))
        d = X[:, :, None] - self.m.T
        var = np.exp(
            - 0.5 * self.ndim / self.beta
            - 0.5 * self.nu * np.einsum('nik, kij, njk -> nk', d, self.W, d)
        )
        r = pi*Lambda*var
        r /= r.sum(axis=-1, keepdims=True)
        r[np.isnan(r)] = 1e-9
        self.pi = pi
        self.Lambda = Lambda
        self.r = r
        return r
    
    def is_converged(self):
        if self._param is None:
            self._param = [self.m, self.W.ravel(), self.alpha, self.beta, self.nu]
            return False
        param = [self.m, self.W.ravel(), self.alpha, self.beta, self.nu]
        if np.array([np.allclose(self._param[i], 
                                 param[i],
                                 rtol=1e-3) for i in range(len(param))]).all():
            return True
        else:
            self._param = param
            return False
        
        
    def maximization(self, X, r):
        Nk = np.sum(r, axis=0)
        xk = 1/Nk * X.T.dot(r)
        dev = X[:,:,None] - xk
        Sk = 1/Nk * np.einsum("nik,njk->ijk", dev,r[:,None,:]*dev)
        
        self.beta = self.beta0 + Nk
        self.alpha = self.alpha0 + Nk
        self.m = ((self.beta0 * self.m0[:, None] + Nk * xk) / self.beta).T
        W0inv = np.tile(np.linalg.inv(self.W0), (self.num_clusters, 1, 1))
        d = (xk - self.m0[:,None])
        var = np.einsum("ik,jk->kij",d,d)
        bn = self.beta0*Nk/(self.beta0+Nk)
        self.W = W0inv + (Nk.reshape((1,1,-1))*Sk).transpose((2,0,1))+\
                     (bn.reshape((-1,1,1))*var)
        self.W = np.linalg.inv(self.W)
        self.nu = self.nu0 + Nk
        self.Sk = Sk.transpose((2,0,1))
        self.xk = xk
        self.Nk = Nk
        
    def elbo(self):
        pass
    
    def fit_init(self, X):
        self.elbos = []
        self.num_samples, self.ndim = X.shape
        
        self.m0 = np.zeros(self.ndim)
        self.W0 = np.eye(self.ndim)
        self.beta0 = 1.
        self.alpha0 = 1.
        self.nu0 = self.ndim
        
        self.m = X[np.random.choice(self.num_samples, self.num_clusters,replace=False)]
        self.W = np.tile(np.eye(self.ndim), (self.num_clusters, 1, 1))
        self.alpha = np.ones((self.num_clusters))
        self.beta = np.ones((self.num_clusters)) * self.beta0
        self.nu = np.ones(self.num_clusters) * self.ndim
        
    def fit_1(self, X, i):
        r = self.expectation(X)
        self.r = r
        self.maximization(X, r)
        # self.elbos.append(self.elbo())
        if self.is_converged():
            print('Iter:{} all parameter converged'.format(i+1))
            return True
        return False

    def fit(self, X, num_iter=1000):
        self.fit_init(X)
        for i in range(num_iter):
            if self.fit_1(X,i):
                break
        else:
            print("parameters may not be converged yet")
                
    def predict(self, X):
        s_alpha = self.alpha.sum()
        L = np.array([((self.nu + 1 -self.ndim) * self.beta / (1 + self.beta))[k] *\
                      self.W[k]for k in range(len(self.W))])
        St =lambda k : multivariate_student_t(X,self.m[k],
                                              np.linalg.inv(L[k]), 
                                              self.nu[k] + 1 - self.ndim + 1)
        prob = np.sum(np.array([self.alpha[k] * St(k)
                                for k in range(len(self.W))]), axis=0)/s_alpha
        return prob
    
    def classify(self, X):
        r = self.expectation(X)
        return np.argmax(r, axis=1)
    

テンソルを扱う際のアインシュタインの縮約記法が非常に便利であることが見て取れるのではないかと思います.

実際にクラスタリングを行ってみましょう.

def plot(model, X,title='Variational Gaussian Mexture Model',ndiv=100,save=None):
    plt.figure(figsize=(7,7))
    X0_range = X[:,0].min()-1, X[:,0].max()+1
    X1_range = X[:,1].min()-1, X[:,1].max()+1
    Y = model.classify(X)

    x, y = np.meshgrid(np.linspace(*X0_range,num=ndiv),
                       np.linspace(*X1_range,num=ndiv))
    plt.scatter(x=X[:, 0], y=X[:,1], alpha=0.6,c=cm.jet(Y / Y.max()))
    Z = model.predict(np.array([x.ravel(), y.ravel()]).T).reshape(x.shape)
    CS = plt.contour(x, y, Z, cmap='jet', alpha=0.7)
    plt.clabel(CS, inline=1, fontsize=10)
    for mu,var in zip(model.m, model.W):    
        plt.scatter(mu[0], mu[1],
                    marker='o',
                    facecolors='red',
                    edgecolors='black',
                    s=70, alpha=0.7)

    plt.title(title)
    if save is not None:
        plt.savefig(save)
plt.figure(figsize=(7,7))
X, y = datasets.make_blobs(n_features=2, centers=4, n_samples=300)
plt.scatter(X[:,0],X[:,1])

model = VariationalMixtureGaussian(7)
model.fit(X, 10000)
plot(model, X, ndiv=50)

収束の様子をアニメーションにしたものは以下のようになる.


コーヒーブレイク

markdownエディタとしておすすめなのがTyporaです.

いままで触ってきた中でかなり使いやすいエディタだと思います.

qiita.com

www.doz3.info

この記事がとても詳しく,プレビュー表示画面というものが存在しなく,表やTODOの追加などを視覚的に行えて(Dropbox Notepaperがエディタとしては近い)

CSSでスタイルを変更でき,論文チックな見やすい見た目に変更することができます.

下書きをいままでHackMDで書いていたのですが操作性可視性ともにこちらに軍配が上がっており,これからはこれをメインに使って行こうかと思います.


はかどるよ

映像のための音楽~平沢進サウンドトラックの世界

パプリカまだみてないのでいつか見てみたい

🐜本の進捗(ベルマンフォード/ダイクストラ/クラスカル)::9/14

得るものが大きい.

2-5 グラフ

頂点と辺からなるデータ構造.頂点と辺の集合を$V,E$であるようなグラフを$G=(V,E)$と表し,2点$u,v$を結ぶ辺を$e=(u,v)$のように表す.

グラフには結びつきに方向性のある有向グラフと,方向性のない無向グラフがある.

木構造は閉路を持たない無向グラフであり,かつ連結であるものを指す.

代わりに連結でないものは森と呼ばれる.

有向グラフに対し,とある頂点$v$から出ていく辺の集合と,入ってくる辺の集合はそれぞれ,$\delta_+(v),\delta_-(v)$とあらわされ,$|\delta_+(v)|,|\delta_-(v)|$をそれぞれ,頂点$v$の出次数,入次数という.

閉路を持たない有向グラフはDAG(Directed Acyclic Graph:有向非巡回グラフ)と呼ばれる.

実装時の方法

行列を使って辺の情報を記憶する方法と,各々の点にリストを与え,辺を記憶する方法があるが,前者はメモリを,疎グラフであっても常に頂点数の二乗だけ必要となってしまうという欠点がある.

実装

C++の勉強兼復習を兼ねて,classでグラフを定義するとどうなるかを考える.

クラス設計等の知識が皆無なので,この方法で良いのかはわからない.

#include <bits/stdc++.h>
using namespace std;
int const nil = (1 << 31);

template <typename T>
class Vertex{
    public:
      int param=nil;
      vector< pair<Vertex *, T> > e;
      int size(){
          return e.size();
      }
};

template <typename T>
class Graph{
    public:
    vector<Vertex<T> *> V;
    Vertex<T>* getVertex(int idx){
        return V[idx];
    }
    int getColor(int idx){
        return V[idx]->param;
    }
    Graph(int v_size){
        for (int i = 0; i < v_size;i++){
            V.push_back(new Vertex<T>());
        }
    }
    int size(){
        return V.size();
    }
    void colorize(int idx, int color){
        V[idx]->param = color;
    }
    void colorize(Vertex<T>* vertex, int color){
        vertex->param = color;
    }
    bool isColorize(int idx, int color){
        if(V[idx]->param==nil)
            return false;
        V[idx]->param = color;
        return true;
    }
    void unite(int v_from, int v_to,T weight=1, bool digraph=true){
        V[v_from]->e.push_back(make_pair(V[v_to], weight));
        if(!digraph){
            V[v_to]->e.push_back(make_pair(V[v_from], weight));
        }
    }
};

彩色問題に対応するため,頂点を簡単なクラスとしてつながりと色を保持するようにする.グラフ自体をクラスとして,結合や色塗りが関数で行えるようにする.

彩色問題の探索の例として以下の問題を考えてみる.

頂点数$N$の無向グラフのすべての頂点に対し,隣り合う頂点同士が違う色になるように2色で色を塗ることができるかを出力せよ.

入力が辺の数を$K$として,

$N~K\\$

$v_1~u_1\\$

$...\\$

$v_k~u_k$

と与えられると考えると,深さ優先探索を用いて,以下のように解ける.

bool dfs(Graph<int>* g, Vertex<int>* v, int c){
    g->colorize(v, c);
    auto edges = v->e;
    for (auto u : edges)
    {
        if(v->param == u.first->param)
            return false;
        if(u.first->param==nil && !dfs(g, u.first, -c))
            return false;
    }
    return true;
}

int main(){
    int N, E;
    cin >> N >> E;
    Graph<int> *g = new Graph<int>(N);
    for (int i = 0; i < E; i++){
        int a, b;
        cin >> a >> b;
        g->unite(a, b, 1, false);
    }
    for(auto v : g->V){
        if(v->param==nil){
            if(!dfs(g, v, 1)){
                cout << "No" << endl;
                return 0;
            }
        }
    }
    cout << "Yes" << endl;
    return 0;
}

最短路問題

Bellman-Ford法

重み付き有向グラフにおける単一始点の最短経路問題を解くラベル修正アルゴリズム[1]の一種である。各辺の重みは負数でもよい。辺の重みが非負数ならば優先度付きキューを併用したダイクストラ法の方が速いので、ベルマンフォード法は辺の重みに負数が存在する場合に主に使われる。

(wikipediaより)

全辺を走査的に見ていって,始点からの最短距離を$d[i]$としたとき, $d[i]=\min\{d[j]+cost(j\rightarrow i )|e=(j,i)\in E\}$が最短経路であることを利用して,更新を行う.負の閉路が存在しないと,高々更新は$|V|-1$回で済むので,全体としては$O(|V|\cdot|E|)$の計算量となる.

Dijkstra法

グラフ理論における辺の重みが非負数の場合の単一始点最短経路問題を解くための最良優先探索によるアルゴリズムである。辺の重みに負数を含む場合はベルマン-フォード法などが使える。辺の重みが全て同一の非負数の場合は幅優先探索が速く、線形時間で最短路を計算可能である。

(wikipediaより)

多分かなり有名なアルゴリズムなので説明はいらないよね.

  1. 各地点までの距離を未確定とし、とりあえず無限大としておきます。
  2. 始点の距離を0とおきます。
  3. 未確定の地点の中から、距離が最も小さい地点を選んで、 その距離を「その地点までの最短距離」として確定します。
  4. 今確定した地点から「直接つながっている」かつ 「未確定である」地点に対して、今確定した場所を経由した場合の 距離を計算し、今までの距離よりも小さければ書き直します。
  5. 全ての地点が確定すれば終了です。そうでなければ3へ。

優先度付きキューを用いることで,最小値の取り出しが容易に行なえ,計算時間の短縮($O(|V|^{2})\rightarrow O(|E|\log|V|)$)が見込めます.

最小全域木

まず,全域木(Spanning tree)とは連結グラフの全ての頂点とそのグラフを構成する辺の一部分のみで構成される木のこと.つまり,連結グラフから適当な辺を取り除いていき,閉路をもたない木の形にしたものが全域木となる.ここで,グラフの各辺に重みがある場合,重みの総和が最小になるように辺を選んで作った全域木のことを最小全域木(Minimum spanning tree)という.

(こちらのサイトより)

最小全域木を求めるには主に以下の2つのアルゴリズムがある.

  1. Kruskal法
  2. Prim法

うち,Kruskal法について実装した.

Kruskal法

まず,頂点全体をそれぞれ1つの木とみなす.(木の数=頂点の数の森を作る)

全辺の中から,重みが最小のものを持ってきて,全体のコストに付け足す.

全体が一つの木になるまでこれを繰り返す.ただし,辺が結ぶ2点が同一の木の要素であるならば,閉路となるので飛ばす.

2を結ぶ...

3を結ぶ...

4を結...ばない(頂点1と4が同一の木構造に属するため)

5を結んで全体が1つの木になる.全体のコストの和は13となる.

このアルゴリズムを実装するに当たって肝になるのが,頂点が同一の木に属するかの判定である.これは,UnionFindを用いるのが非常に相性が良い.UnionFindの計算量は対数より小さいアッカーマンの逆関数$O(\alpha(|V|))$となるので全体の計算量は,ソートアルゴリズムの計算量で決まる.

実装

上記3つのアルゴリズムを実装した.

入力は最小全域木の説明で用いたものと同一とし,経路探索の場合は頂点0から5への最短距離(=7)の出力を目的とした.

#include <bits/stdc++.h>
using namespace std;
static const int nil = (1 << 31);
static const int inf = (1 << 30);

class UnionFind{
    private:
    vector<int> arr;
    public:
    UnionFind(int size): arr(size, -1){}
    int root(int node){
        return arr[node] < 0 ? node : arr[node] = root(arr[node]);
    }
    int size(int x){
        return -arr[root(x)];
    }
    bool unionSet(int x, int y){
        x = root(x);
        y = root(y);
        if(x!=y){
            if(arr[x] < arr[y])
                swap(x, y);
            arr[x] += arr[y];
            arr[y] = x;
        }
        return x != y;
    }
    bool findSet(int x, int y){
        return root(x) == root(y);
    }
};

template <typename T>
class Vertex{
    public:
      int idx;
      int param = nil;
      int dist = inf;
      vector<pair<Vertex *, T> > e;
      int size(){
          return e.size();
      }

      bool operator ==(const Vertex* o){
          return dist == o->dist;
      }
      bool operator !=(const Vertex* o){
          return dist != o->dist;
      }
      bool operator >(const Vertex* o){
          return dist > o->dist;
      }
      bool operator >=(const Vertex* o){
          return dist >= o->dist;
      }
      bool operator <(const Vertex* o){
          return dist < o->dist;
      }
      bool operator <=(const Vertex* o){
          return dist <= o->dist;
      }
};

template <typename T>
class Graph{
    public:
    vector<Vertex<T> *> V;
    typedef tuple<Vertex<T> *, Vertex<T> *, T> Edge;
    vector<Edge> E;
    Graph(int v_size){
        for (int i = 0; i < v_size;i++){
            auto v = new Vertex<T>();
            v->idx = i;
            V.push_back(v);
        }
    }
    static bool comp(const Edge &e1, const Edge &e2)
    {
        return get<2>(e1) < get<2>(e2);
    }
    Vertex<T> *getVertex(int idx)
    {
        return V[idx];
    }
    int getColor(int idx){
        return V[idx]->param;
    }
    int size(){
        return V.size();
    }
    void colorize(int idx, int color){
        V[idx]->param = color;
    }
    void colorize(Vertex<T>* vertex, int color){
        vertex->param = color;
    }
    bool isColorize(int idx, int color){
        if(V[idx]->param==nil)
            return false;
        V[idx]->param = color;
        return true;
    }
    void unite(int v_from, int v_to,T weight=1, bool digraph=true){
        E.push_back(make_tuple(V[v_from], V[v_to], weight));
        V[v_from]->e.push_back(make_pair(V[v_to], weight));
        if(!digraph){
            E.push_back(make_tuple(V[v_to], V[v_from], weight));
            V[v_to]->e.push_back(make_pair(V[v_from], weight));
        }
    }
    void dijkstra(int idx){
        priority_queue<Vertex<T> *, vector<Vertex<T> *>, greater<Vertex<T> *> > pq;
        Vertex<T>* s = getVertex(idx);
        s->dist = 0;
        pq.push(s);
        while (!pq.empty())
        {
            Vertex<T>* v = pq.top();
            pq.pop();
            for (auto u : v->e){
                if(u.first->dist > v->dist + u.second){
                    u.first->dist = v->dist + u.second;
                    pq.push(u.first);
                }
            }
        }
    }
    void bellmanford(int idx){
        Vertex<T>* s = getVertex(idx);
        s->dist = 0;
        while(true){
            bool update = false;
            for(auto e : E){
                Vertex<T> *from = get<0>(e);
                Vertex<T> *to = get<1>(e);
                T cost = get<2>(e);
                if(from->dist != inf && to->dist > from->dist + cost){
                    to->dist = from->dist + cost;
                    update = true;
                }
            }
            if(!update)
                break;
        }
    }
    T kruskal(){
        sort(E.begin(), E.end(), Graph::comp);
        UnionFind uf = UnionFind(size());
        int res = 0;
        for (auto e : E){
            auto u = get<0>(e);
            auto v = get<1>(e);
            if(!uf.findSet(u->idx, v->idx)){
                uf.unionSet(u->idx, v->idx);
                res += get<2>(e);
            }
        }
        return res;
    }
    void stateInit(){
        for(auto v : V){
            v->dist = inf;
            v->param = nil;
        }
    }
};


int main(){
    int N, E;
    cin >> N >> E;
    Graph<int> *g = new Graph<int>(N);
    for (int i = 0; i < E; i++){
        int a, b, c;
        cin >> a >> b >> c;
        g->unite(a, b, c, false);
    }
    g->dijkstra(0);
    cout << g->getVertex(N - 1)->dist << endl;
    g->stateInit();
    g->bellmanford(0);
    cout << g->getVertex(N - 1)->dist << endl;
    g->stateInit();
    int v = g->kruskal();
    cout << v << endl;

    return 0;
}

実行例

入力:
6 9
0 1 2
1 3 3
3 5 7
0 2 5
2 4 9
4 5 1
1 4 4
1 5 6
3 4 2

出力:
7
7
13

🐜本の進捗(優先度付きキューとUnionFind木)::9/13

いままでまともにアルゴリズムの勉強をしてこなかったので得るものが大きいです.

2-4 データ構造

二分木

すべてのノードの子の数が2以下の木

プライオリティキュー

数の追加と最小値の取り出しが行えるデータ構造を指す

二分木ヒープ

プライオリティキューを二分木を用いて実装したもの

子の数字>親の数字を満たす二分木で更新と削除のたびに左記条件を満たすまで交換を繰り返す.

実装方法

#include <bits/stdc++.h>
using namespace std;
#define REP(i, n) for (int i = 0; i < n; i++)

template <typename T>
class Heap
{
  public:
    int sz = 0;
    T array[100001];
    Heap()
    {
    }
    void push(T x)
    {
        int i = sz++;
        while (i > 0)
        {
            int p = (i - 1) / 2;
            if (array[p] <= x)
                break;
            array[i] = array[p];
            i = p;
        }
        array[i] = x;
    }

    T pop()
    {
        T ret = array[0];

        T x = array[--sz];
        int i = 0;
        while (i * 2 + 1 < sz)
        {
            int a = i * 2 + 1, b = i * 2 + 2;
            if (b < sz && array[b] < array[a])
                a = b;
            if (array[a] >= x)
                break;
            array[i] = array[a];
            i = a;
        }
        array[i] = x;
        return ret;
    }
};

int main(void)
{
    Heap<float> heap;
    heap.push(1.6);
    heap.push(7.9);
    heap.push(-44.1);
    heap.push(842.0);
    heap.push(0.0);
    REP(i, 5)
    {
        cout << heap.pop() << endl;
        /*
            -44.1
            0
            1.6
            7.9
            842
        */
    }
    return 0;
}

二分木なので親ノードのインデックスは(i-2)/2で取得でき,子ノードはi*2+1, i*2+2で取得が可能であることに注意しよう!!!!

$\text{function}~~push(Num: x) \rightarrow Nil$

  1. $arr = Array.new();sz = 0;//initialize$
  2. $i = sz++; // i \leftarrow \text{Self Node Index}$
  3. $\text{while}(i > 0):$
    1. $p = (i-1)/2; // p \leftarrow \text{Parent Node Index}$
    2. $\text{if}(arr[p] \leq x) \text{~then~~break;}$//逆転していたら抜ける
    3. $arr[i]=arr[p];~i=p;$
  4. $arr[i]=x$

$\text{function}~pop() \rightarrow Num$

  1. $retval = arr[0];$//返り値(最小値)
  2. $x = arr[-1];$//根に持ってくる値
  3. $i=0;$
  4. $\text{while}(i*2+1<sz):$//子ノードが存在する
    1. $l=i\times2+1,r=i\times2+2;$
    2. $\text{if}(r<sz$ $and$ $arr[b]<arr[a])\text{then}$ $a\leftarrow b;$
    3. $\text{if}(arr[a] \geq x)\text{thenbreak;}$
    4. $arr[i]=arr[a];~i=a;$
  5. $arr[i]=x;$
  6. $\text{return}~~retval;$

c++ のSTLにもプライオリティキューはある.そのため,イチから自分で実装するという機会はほぼないが,STLの方は最大値から取り出されるので注意が必要.

DS-1(POJ 2431)

壊れた燃料タンクを直すために距離が $L$ だけ離れている街へ向かう. タンクには最初 $P$ だけ燃料があるが,単位距離進む毎に1単位分の燃料を失う. 街までの道には $N$ 個の燃料補給所がある. $i$ 番目の補給所は故障地点から $A_i$ の位置にあり,そこで $B_i$ の燃料を補給することができる. タンクに保存できる燃料の量に上限は無い. 街まで進むために必要な最小の燃料補給の回数を答えよ. ただし,どのように進んでも街まで辿り着けない場合は $-1$ と出力せよ.

(POJ本家は街からの距離を$A_i$としていたが,蟻本では簡単のため故障地点からの位置に変えてある.)

ガソリンスタンドに寄って補給するのではなく,ガソリンスタンドに寄ったらそのガソリンスタンドで燃料を補給する権利を手に入れると考える.

つまり,通過したガソリンスタンドの補給量をPQに記憶し,燃料が尽きた地点で,立ち寄ったガソリンスタンドの中で最も補給量が大きいものを選べば良いことがわかる.

#include <queue>
#include <iostream>

using namespace std;
#define REP(i, n) for (int i = 0; i < n; i++)

int main()
{
    int N, L, P;
    int A[10001], B[10001];
    cin >> N;
    REP(i, N)
    {
        cin >> A[i]] >> B[i];
    }
    cin >> L >> P;

    A[N] = L;
    B[N] = 0;
    N++;
    priority_queue<int> pq;

    int ans = 0, pos = 0, tank = P;
    REP(i, N)
    {
        int d = A[i] - pos;
        while (tank - d < 0)
        {
            if (pq.empty())
            {
                cout << "-1" << endl;
                return 0;
            }
            tank += pq.top();
            pq.pop();
            ans++;
        }
        tank -= d;
        pos = A[i];
        pq.push(B[i]);
    }
    cout << ans << endl;
    return 0;
}

Union-Find木

グループ分けを管理するデータ構造の一種.次のことが効率的に行える.

  • 要素$a$と要素$b$が同じグループに属するかを調べる
  • 要素$a$と要素$b$のグループを併合する.分割は出来ない

一つの木を一つのグループとして扱う.

併合

木同士の根に辺を張ることで,グループを併合することができる.

要素判定

根が同じ要素同士は同一グループとみなすことができる.

実装上の注意

木構造を扱うので,深さに偏りが発生すると計算量が増大してしまう.そこで,各木の深さを記憶してとき,深さの浅い木を深い方に辺を張ることで偏りの少ない木とすることができる.

また,各ノードが根を辿った場合,辺を直接根に張り直すことで根からの深さを抑えることができる.

実装例

http://www.prefield.com/algorithm/container/union_find.html

class UnionFind{
    private:
    vector<int> arr;
    public:
    UnionFind(int size): arr(size, -1){}
    int root(int node){
        return arr[node] < 0 ? node : arr[node] = root(arr[node]);
    }
    int size(int x){
        return -arr[root(x)];
    }
    bool unionSet(int x, int y){
        x = root(x);
        y = root(y);
        if(x!=y){
            if(arr[x] < arr[y])
                swap(x, y);
            arr[x] += arr[y];
            arr[y] = x;
        }
        return x != y;
    }
    bool findSet(int x, int y){
        return root(x) == root(y);
    }
};
DS-2(POJ 1182)

$N$匹の動物に対し$1..N$の番号が振られており,各々の動物は3つの$A,B,C$に 分けられる.このグルーピングには以下のような強弱関係がある.$A\rightarrow B \rightarrow C \rightarrow A$.

次の2種類の情報が$k$個与えられる.

  1. $x$と$y$は同じ種類である
  2. $x$は$y$を捕食する

入力を逐次処理した際,矛盾する記述は何個あるか.ただし矛盾した記述は以降捨てると考える.

int main(void)
{
    bool debug = true;
    int N, k;
    int T[50001], X[50001], Y[50001];
    cin >> N >> k;
    REP(i, k)
    {
        scanf("%d %d %d", &T[i], &X[i], &Y[i]);
    }
    UnionFind uf = UnionFind(3 * N);

    int ans = 0;
    REP(i, k){
        int t = T[i];
        int x = X[i] - 1;
        int y = Y[i] - 1;
        if (x < 0 || x >= N || y < 0 || y>=N){
            ans++;
            if (debug)
                cout << "[*] Invalid Data (Out of length): x = " << x << ", y = " << y << endl;
            continue;
        }

        // Same Group
        if(t==1){
            if(uf.findSet(x, y+N) || uf.findSet(x, y+2*N)){
                ans++;
                if (debug){
                    cout << "[*] Invalid 'Same kind' input: x = ";
                    cout << x << ", y = " << y << endl;
                }
            }else{
                uf.unionSet(x, y);
                uf.unionSet(x + N, y + N);
                uf.unionSet(x + 2 * N, y + 2 * N);
                if (debug){
                    cout << "[*] Valid 'Same kind' input: x = ";
                    cout << x << ", y = " << y << endl;
                }
            }
        }else if(t == 2){ // Eat Group
            if(uf.findSet(x, y+N) || uf.findSet(x, y+2*N)){
                ans++;
                if (debug){
                    cout << "[*] Invalid 'Eat kind' input: x = ";
                    cout << x << ", y = " << y << endl;
                }
            }else{
                uf.unionSet(x, y + N);
                uf.unionSet(x + N, y + 2*N);
                uf.unionSet(x + 2 * N, y);
                if (debug){
                    cout << "[*] Valid 'Eat kind' input: x = ";
                    cout << x << ", y = " << y << endl;
                }
            }
        }else{
            ans++;
            cout << "[*] Invalid Type: t = " << t;
        }
    }
    cout << ans << endl;
    return 0;
}

3*Nのノードを持つUnionFindを定義しておいて,$0..N-1$をグループA,$N..2N-1$をグループB,$2N..3N-1$をグループCとすることで同族ならば実線のように,捕食関係にあるならば破線のように互いを結び付けることで,各々の関係を記録することができる.

初心に帰って変分ベイズ(1)::8/31

1. 確率モデルの近似法

確率モデルでは,観測データ${\bf X}$に対するモデルの潜在変数${\bf Z}$の事後分布$p({\bf Z}|{\bf X})$を求めることが目的です.しかしながら,「確率」として取り扱うには総和が1である(=正規化する)必要があります.ベイズの定理においては分母に相当すると言えて, $$ p({\bf Z}|{\bf X})=\frac{p({\bf X}|{\bf Z})p({\bf Z})}{p({\bf X})} $$ というベイズの式から,分母の周辺尤度$p({\bf X})$は分子を$d{\bf Z}$で積分(=周辺化)を行うことで得られて $$ p({\bf X})=\int p({\bf X}|{\bf Z}) p({\bf Z})d{\bf Z} $$ という式になります.

しかしながら,尤度関数($p({\bf X}|{\bf Z})$)事前分布($p({\bf Z})$)がどちらも正規分布などという簡単な条件下では解析的に解けるものの,複雑な一般の問題ではこの計算を行えない場合が多いです. というのも,潜在変数空間の次元数が膨大であったり,積分が閉じた形式*1の解を持たないといったことが主な理由です.

解析的に求まらない問題に対しては何を行うかと言えば,やはり近似です.

よく使われる近似法としては

  • マルコフ連鎖モンテカルロ(MCMC)
  • 変分ベイズ法(変分推定法)

があり,前者が確率的近似に対して後者が決定的近似です.

MCMCについての詳細は後に述べるとしてそれぞれ以下のような長所と短所が存在します.

長所 短所
MCMC $N\rightarrow\infty$ならば厳密解が得られる 計算量が嵩む
変分ベイズ法 計算が少ない 厳密解は得られない

変分ベイズ法は複雑な事後分布を分解し近似を行うというスタンスで学習を行います.

その理論の導出の前に,導出の際に用いるKLダイバージェンスについて今だ理解が足りてなかったのでまとめます*2

2. エントロピーとKLダイバージェンス

2.1. エントロピー

情報学における(離散分布の)エントロピーは以下のように定義でき, $$ I(X)=-\sum_{i}p(i)\log{p(i)} $$ 名の示す通り,統計力学のエントロピーと同様,確率変数$X$の取りうる「乱雑さ」を定量的に表す指標といえます.

例えば,6面体サイコロ$X=[1,2,3,4,5,6]$を考えたとき出目の確率が

  • $p(X)=[1/5, 1/10, 1/10, 1/5, 1/5, 1/5]$
  • $q(X)=[1/2,1/4,1/8,1/16,1/32,1/32]$
  • $r(X)=[1/6,1/6,1/6,1/6,1/6,1/6]$

であるします.

エントロピーはこの場合出目が偏らないほど予測が行いづらいので,いわば出目の予測のしにくさを表す量ともいえます.

それぞれのエントロピーをpythonで求めてみます.

import numpy as np
import matplotlib.pyplot as plt

x = np.arange(1,7)
q = np.array([1/2, 1/4, 1/8, 1/16, 1/32, 1/32])
p = np.array([1/5, 1/10, 1/10, 1/5,  1/5,  1/5])
r = np.ones(6)/6

# 2つの確率分布を可視化する関数です
def plot(p, q):
    m, s, b = plt.stem(x, p)
    plt.setp(m, 'color', 'g')
    plt.setp(s, 'color', 'g')
    m, s, b = plt.stem(x, q)
    plt.setp(s, 'linestyle', '-.', 'alpha', 0.6)
    plt.show()

# エントロピーの定義
H = lambda p: -np.sum(p*np.log(p))

print(H(p)) # 1.74806...
print(H(q)) # 1.3429...
print(H(r)) # 1.79175...
pとrの分布

qとrの分布

予測通り同様に確からしいサイコロのエントロピーが最も大きいことが分かりました.

また,pとrをみると定義式から見ても明らかですが似た分布ではエントロピーの値が近いことも言えます.

2.2. KLダイバージェンス

二つの確率分布$p,q$に対して距離的な解釈を与えたいとします.すなわち,本来の分布$p$から,実際の分布$q$がどの程度離れているかを考えます.

そこで,下式より定義されるKLダイバージェンスを導入します.

$$ D_{KL}(p||q)=\sum_{i}p_i\log{\frac{p_i}{q_i}} $$

これは変形すると$-p_{i}\log{q_{i}}$と$-p_{i}\log{p_{i}}$の差に分解されるので「$q$に従ってデータが発生したと想定して計算したエントロピー」と「実際に$p$に従って観測されたデータによるエントロピー」の差であると解釈できます*3

また,定義からもわかるように可換ではありません.

サイコロの例で実際に計算してみると,

D_KL = lambda p, q:np.sum(p*np.log(p/q))
print(D_KL(p, r)) # 0.043692...
print(D_KL(r, p)) # 0.048727...
print(D_KL(q, r)) # 0.448786...

分布が似ている$p,r$では低い値を取ることが確認できます.

もちろん,KLダイバージェンスは総和を適切に積分に置き換えることで連続な確率変数に対しても適用することが可能です.

このグラフは$p$が標準正規分布に従い,$q(x)$を$q(x)\sim\mathcal N(x;\mu_{q},\sigma_{q})$として横軸を$\mu_{q}$,縦軸を$\sigma_{q}$としたときのKLダイバージェンスをカラーマップで表したものです.$\mu_{q}=0,\sigma_{q}=1$の際に最も青く表示されていることが見て取れます.

定義からもわかるようにKLダイバージェンスは非負で0となるのは分布が完全に等しいときのみという性質があります.

3. 変分ベイズ

3.1. 変分下界

KLダイバージェンスが分布の類似度を測るものだということが分かりました.なので,確率モデルの事後分布$p({\bf Z}|{\bf X})$を分布$q({\bf Z})$で近似させたいとき, $$ D_{KL}(q({\bf Z})||p({\bf Z}|{\bf X})) $$ を最小化させることでその目的が達成されることが容易にわかるかと思います.

それを踏まえたうえで,モデルの対数周辺尤度$\ln p({\bf X})$(対数周辺尤度はしばしエビデンスと呼ばれます)を考えます.

また,この際一般の期待値に関しイェンセンの不等式が成り立つことを利用します

イェンセンの不等式 $$ f(\mathbb E[X])\leq \mathbb E[f(X)] $$

これを利用して

$$ \begin{eqnarray} \ln p({\bf X})&=&\ln\int p({\bf X},{\bf Z})d{\bf Z}\\ &=&\ln\int p({\bf X},{\bf Z})\frac{q({\bf Z})}{q({\bf Z})}d{\bf Z}\\ &=&\ln\mathbb E_{q}\left[\frac{p({\bf X},{\bf Z})}{q({\bf Z})}\right]\\ &\geq&\mathbb E_{q}\left[\ln\frac{p({\bf X},{\bf Z})}{q({\bf Z})}\right]\\ &=& \int q({\bf Z})\ln\frac{p({\bf X},{\bf Z})}{q({\bf Z})}d{\bf Z} \end{eqnarray} $$

この下界を$\mathcal L(q)$のようにして汎関数を定めます.$\mathcal L$は変分下界と呼ばれ,対数周辺尤度がエビデンスと呼ばれることからELBO(Evidence Lower BOund)と呼ばれます.

さて,つぎに対数周辺尤度からこの変分下界を引いてみます.

$$ \begin{eqnarray} \ln p({\bf X})-\mathcal L(q)&=&\int q({\bf Z})d{\bf Z}\ln p({\bf X})-\int q({\bf Z})\ln\frac{p({\bf X},{\bf Z})}{q({\bf Z})}d{\bf Z}\\ &=&\int q({\bf Z})\ln \frac{q({\bf Z})}{p({\bf Z}|{\bf X})}d{\bf Z}\\ &=&D_{KL}(q({\bf Z})||p({\bf Z}|{\bf X})) \end{eqnarray} $$

とKLダイバージェンスが導き出せます.周辺対数尤度は定数であるので,変分下界の最大化とKLダイバージェンスの最小化は等価であることが分かります.

3.2. 分解による近似分布の計算

上記の式変形で汎関数$\mathcal L(q)$が得られました.ここで,汎関数の極値を求めるにあたってオイラーの方程式*4を使います.


オイラーの方程式

汎関数 $$ I[y]\equiv\int \frac{f(x,y,y')}{q(y)}dx $$

の停留値を求めるには, $$ \frac{\partial f}{\partial y}-\frac{d}{dx}\left(\frac{\partial f}{\partial y'}\right)=0 $$

を解けばよい.


この場合は$y'$が$f$に含まれていないので単に, $$ \frac{\partial f(x,y)}{\partial y(x)}=0 $$

を解けばよいとといえます.

さて,このオイラーの方程式をこのモデルに用いて汎関数である変分下界の極値を求めます.

$q({\bf Z})$を排反な分布に分解して${\bf Z_i}(i=1,2,...,M)$と書くとします.このような分解された形を用いて近似を行うことは物理学の文脈では平均場近似(Mean-Field Approximation)と呼ばれ,特に事後分布をこの方法で近似する方法は変分ベイズ法(Variational Bayes)と呼ばれます.

$$ q({\bf Z})=\prod_{i=1}^{M}q_i({\bf Z_i}) $$

この下では変分下界は $$ \mathcal L(q)=\int \prod_{i=1}^{M}q_{i}({\bf Z_i})\left(\ln p({\bf X},{\bf Z})-\sum_{i=1}^{M}\ln q_{i}({\bf Z}_i)\right)d{\bf Z} $$ と書けます.簡単のために2つの場合を考えることにすると, $$ \mathcal L(q)=\iint q({\bf Z}_1)q({\bf Z}_2)\ln\frac{p({\bf X},{\bf Z})}{q({\bf Z}_1)q({\bf Z}_2)} d{\bf Z}_1d{\bf Z}_2 $$ となります.ここで,$q_2({\bf Z}_2)$を固定して考えると$q_1({\bf Z}_1)$の変分下界は,

$$ \mathcal L(q_1)=\int\left( \int q({\bf Z}_1)q({\bf Z}_2)\ln\frac{p({\bf X},{\bf Z})}{q({\bf Z}_1)q({\bf Z}_2)}d{\bf Z}_1\right)d{\bf Z}_2 $$

と式変形ができることから, $$ f(q_1,{\bf Z}_1)=\int q({\bf Z}_1)q({\bf Z}_2)\ln\frac{p({\bf X},{\bf Z})}{q({\bf Z}_1)q({\bf Z}_2)}d{\bf Z}_1 $$ とすると,変分法が適用でき, $$ \frac{\partial f(q_1,{\bf Z}_1)}{\partial q_1({\bf Z}_1)}=0 $$ を解くことで, $$ \ln q_1({\bf Z}_1)=\int q_2({\bf Z}_2)\ln p({\bf X},{\bf Z})d{\bf Z}_2+\text{const} $$ の時に極値を取ることが分かります. これは $$ \ln q_1^{*}({\bf Z}_1)=\mathbb E_{q_2}[\ln p({\bf X},{\bf Z})]+\text{const} $$ とも解釈できます.

$q_2$も同様, $$ \ln q_2^{*}({\bf Z}_2)=\mathbb E_{q_1}[\ln p({\bf X},{\bf Z})]+\text{const} $$

となります.

さらに,一般の$M$分割に対しても同様のことが言え,同時分布を自分以外の分解された近似分布で期待値を取ることで求めることができます.すなわち,$\mathbb E_{i\neq j}$を$i\neq j$であるすべての${\bf Z}_i$による分布$q_i$の期待値とみなすと $$ \ln q^{*}_{j}=\mathbb E_{i\neq j}[\ln p({\bf X},{\bf Z})]+\text{const} $$ のように最適解が与えられます.

変分ベイズ法はこの式を基礎として近似を行っていきます.$M$分割をおこなうことで$M$個の連立方程式が出来上がります.しかし,その方程式は他の因子$q_i({\bf Z}_i)(i\neq j)$の積分に依存しているので解析的に難しいため,各因子を適当な値で初期化して因子をそれぞれ,他の因子を固定して改良を重ねていくことで学習を行う,という手法がとられます.

4. 変分混合ガウス分布

(編集後記:非常に長くなりそうなので2回くらいに分ける.)

4.1. 事前分布の定義

変分ベイズが適用できる例としてクラスタリング問題が考えられます.

上のような分布に対して何個のクラスタリングに分けられるかという問題を考えるとき,K-Means法やEMアルゴリズムという決定論的な方法が考えられますが,ベイズ法を適用することで,特定のデータのみを説明する縮退が起こりえないことや,最適な要素(コンポーネント,クラスタ)数$K$を交差検証無しに求められるといった利点があります.

確率的なクラスタリングでは,この問題を次のようにモデル化します.

まず各観測値${\bf x}_n$に対して,対応するコンポーネントをOneHot化したベクトル${\bf z}_n$とします.データ数を$N$,コンポーネント数を$K$,次元数を$D$とすると,${\bf X}=\{{\bf x}_1,{\bf x}_2,...,{\bf x}_N\}$は$N\times D$,${\bf Z}=\{{\bf z}_1,{\bf z}_2,...,{\bf z}_N\}$は$N\times K$の行列となります.

ここで,混合ガウス分布は,コンポーネント数が$K$のとき,以下のように線形重ね合わせで表現することができることから,

$$ \sum_{k=1}^{K}\pi_k\mathcal N({\bf x}|{\bf \mu}_k,{\bf \Sigma}_k) $$

のように$\sum \pi_k=1,0\leq\pi_k\leq1$を満たす係数$\pi_k$でその混合比を表現することができます.

${\bf z}_n$はOneHotベクトル化されていることから,分布$p({\bf z}_n)$は, $$ p({\bf z}_n)=\prod_{k=1}^{K}\pi_k^{z_{nk}} $$ のように書き下すこともできます.したがって,分布$p({\bf Z}|{\bf \pi})$,$(\text{where}~{\bf \pi}:=\{\pi_1,\pi_2,...,\pi_k\})$は, $$ p({\bf Z}|{\bf \pi})=\prod_{n=1}^{N}\prod_{k=1}^{K}\pi_k^{z_k} $$ となります.また, $$ p({\bf x}_n|{\bf z}_n)=\prod_{k=1}^{K}\mathcal N({\bf x}_n|{\bf \mu}_k,{\bf \Sigma}_k)^{z_{nk}} $$ となることから,潜在変数が与えられたもとでの観測データのベクトルの条件付き分布は $$ p({\bf X}|{\bf Z},{\bf \mu},{\bf \Lambda})=\prod_{n=1}^{N}\prod_{k=1}^{K}\mathcal N({\bf x}_n|{\bf \mu}_k,{\bf \Lambda}_k^{-1})^{z_{nk}} $$

となります.ただし,共分散行列ではなく精度行列${\bf \Lambda}$(共分散行列の逆行列)を計算の簡単化のために用いています.

次にベイズ的アプローチを導入するため,パラメータ${\bf \pi},{\bf \mu},{\bf \Lambda}$に事前分布を導入します.まず,$p({\bf \pi})$は多変量のベータ分布に相当するディリクレ分布を仮定します*5

$$ p({\bf \pi})=\frac{\Gamma(\sum \alpha_{0k})}{\prod\Gamma(\alpha_{0k})}\prod_{k=1}^{K}\pi_{k}^{\alpha_{0k}-1} $$ ただし,${\bf \alpha_{0}}$はすべての要素$\{\alpha_{0k}\}$が同じ値のベクトルとした.また上の式の係数はディリクレの正規化定数と呼ばれ,$C({\bf \alpha}_0)$のように表記する.

次に分布$p({\bf \mu},{\bf \Lambda})$は,平均と分散が未知なのでガウスウィシャート分布が適します. $$ p({\bf \mu},{\bf \Lambda})=\prod_{k=1}^{K}\mathcal N({\bf \mu}_k|{\bf m}_0,(\beta_0{\bf \Lambda_k})^{-1})\mathcal W({\bf \Lambda_k}|{\bf W},\nu_0) $$


(以下ポエム)

残りは次回に回すことにします.この例はPRMLでもなかなか長く歯ごたえのある所だと思うので念入りにやっていきたい.いままでVAEなどについても理解が及ばない点が多く,変分下界の存在意義とかそういったものを理論立てて学べる良い教材だった.

次回

fgjiutx.hatenablog.com

*1:閉形式ともいう

*2:タイトル詐欺か?

*3:「想定・実際」は逆で解釈してもよい

*4:[疑問]オイラーさん方程式作りすぎてどれが変分に関わるものなの?

*5:PRML上巻p305を参照するとわかるが,$\pi$が満たすべき条件に完全に合致している

RustでAtcoder Beginner Contest100のD問題を解いた::8/20

D - Patisserie ABC

高橋君はプロのパティシエになり, AtCoder Beginner Contest 100 を記念して, 「ABC洋菓子店」というお店を開いた.

ABC洋菓子店では, $N$ 種類のケーキを売っている.各種類のケーキには「綺麗さ」「おいしさ」「人気度」の $3$つの値を持ち, $i$種類目のケーキの綺麗さ$x_i$は, おいしさは$y_i$, 人気度$z_i$は である.これらの値は$0$以下である可能性もある. りんごさんは, ABC洋菓子店で$M$個のケーキを食べることにした. 彼は次のように, 食べるケーキの組み合わせを選ぶことにした.

  • 同じ種類のケーキを $2$個以上食べない.
  • 上の条件を満たしつつ, (綺麗さの合計の絶対値) + (おいしさの合計の絶対値) + (人気度の合計の絶対値) が最大になるように選ぶ.

このとき, りんごさんが選ぶケーキの (綺麗さの合計の絶対値) + (おいしさの合計の絶対値) + (人気度の合計の絶対値) の最大値を求めなさい.

色々と格闘した結果書けたけどこれでいいのか~という気持ちになってる

こんな書き方でいいのか強い方に見てほしいという気持ちになっている.

D問題

use std::cmp::max;
use std::fmt::Debug;
use std::io::stdin;
use std::str::FromStr;

fn read<T>() -> Vec<T>
where
    T: FromStr,
    <T as FromStr>::Err: Debug,
{
    let mut s = String::new();
    stdin().read_line(&mut s).unwrap();
    s.trim()
        .split_whitespace()
        .map(|c| T::from_str(c).unwrap())
        .collect()
}

#[derive(Clone)]
struct Review {
    x: i64,
    y: i64,
    z: i64,
}

impl Review {
    fn new(r: &[i64]) -> Review {
        Review {
            x: r[0],
            y: r[1],
            z: r[2],
        }
    }

    fn dot(&self, o: &Review) -> i64 {
        o.x * self.x + o.y * self.y + o.z * self.z
    }
}

fn main() {
    let x: Vec<usize> = read();
    let mut select: Vec<Review> = Vec::new();
    for i in 0..8 {
        let v: Vec<i64> = vec![i & 0x4, i & 0x2, i & 0x1]
            .iter()
            .map(|&i| if i == 0 { -1 } else { 1 })
            .collect();
        select.push(Review::new(v.as_slice()));
    }

    let mut reviews: Vec<Review> = Vec::new();
    for _ in 0..x[0] {
        let r: Vec<i64> = read();
        reviews.push(Review::new(r.as_slice()));
    }
    println!(
        "{}",
        select
            .iter()
            .map(|sel| {
                let mut biased_points: Vec<i64> = reviews.iter()
                                                         .map(|v| v.dot(&sel))
                                                         .collect::<Vec<i64>>();
                biased_points.sort_by(|a, b| b.cmp(a));
                biased_points
                    .iter()
                    .take(x[1])
                    .fold(0, |sum, item| sum + item)
                    .abs()
            })
            .fold(0, |sum, item| max(sum, item))
    )
}

絶対値が最大となる場合であるのでアルゴリズム的には各々の評価指標が負に大きいか正に大きいかを全通り探索($2^{3}=8$通り)すればよい.

しかしまあ,Rustの使い方が全く分からずかなり*1苦戦した.

*1:4時間くらい