独立成分分析の実装

はじめに

統計的因果探索を読んでいる. この中で,独立成分分析が出てくるが,そういえばあんまり理解できていなかったので,実装して理解する. 日本語の本も見たけど,結局 Independent Component Analysis: Algorithms and Applicationsがわかりやすかった.

独立成分分析

独立成分分析では,観測 \textbf{x} = A\textbf{s}が与えられたとき,元の確率変数 \textbf{s} = {s_1, s_2, ..., s_N}を復元することを試みる. ここで,独立成分分析では,確率変数  \textbf{s}の要素が独立であるという仮定をおいて,信号を復元する.

以下では,矩形波と制限波からなる \textbf{s}を例にプログラムを書く.

import pylab as plt
%matplotlib inlinee

import numpy as np
from scipy import signal

from sklearn.decomposition import PCA

import torch
from torch.autograd import Variable

np.random.seed(0)

n_samples = 5000
time = np.linspace(0, 8, n_samples)

s1 = np.sin(2 * time)  # Signal 1 : sinusoidal signal
s2 = np.sign(np.sin(3 * time))  # Signal 2 : square signal

S = np.c_[s1,  s2]

S /= S.std(axis=0)  # Standardize data
A = np.array([[1, 1], 
              [0.5, 2],
             ])  # Mixing matrix

X = np.dot(S, A.T)  # Generate observations

元の信号は以下.

plt.plot(S)
plt.title("S")
plt.show()

f:id:ksknw:20180519231649p:plain

この信号に対して,以下のような観測が得られたとする.

plt.plot(X)
plt.title("X")
plt.show()

f:id:ksknw:20180519231705p:plain

独立成分分析では, \textbf{s}の各要素が独立であるという仮定から,復元行列  Wをいい感じに求めて,復元値 \textbf{y} = W\textbf{x}を決める.

中心極限定理と独立

2次元の系列データ \textbf{y} = [y_1, y_2]を復元するとする. このとき, [tex: y_1 = \textbf{w}1^\top A \textbf{s}] ここで, \textbf{w}1は Wの行ベクトル(1行目)である. ここで, \textbf{z}=A\textbf{w} とすると  y_1 = \textbf{z}^\top \textbf{s}であり,復元値  y_1, y_2 は独立な確率変数 \textbf{s}を重み付きで足したものであるといえる.

中心極限定理から独立な確率変数の和は,ガウス分布に近づくことが多い(参考文献では"usually"). このため,復元がうまく行っていない場合には, \textbf{y} \textbf{s}に比べてガウス分布に近づくことが多い. 独立成分分析では,この性質を利用して \textbf{y}がよりガウス分布から離れるように Wを求める.

ガウス分布度合いを測る指標として尖度(kurtosis)がある.尖度はガウス分布に近いときに0となる.  \textbf{y} の分散共分散行列が単位行列のとき,尖度は \mathrm{kurt} ( \textbf{y} ) = \mathbb{ E } [ y^4 ]-3 として求まる. 独立成分分析では,尖度の絶対値を最大化するような Wを勾配法によって求める. ただし, Wの大きさを大きくすると,どこまでの誤差関数を小さくできてしまうので, Wのノルムを1に制約する. 今回は適当に毎ステップでノルムが1になるようにリスケールすることで,制約を満たす解が求まるようにした.

pca = PCA(n_components = X.shape[1], whiten=True)
X = pca.fit_transform(X)
def kurt(Y):
    return torch.mean(Y**4, dim=0) - 3
def kurt_loss(Y):
    # Y の分散共分散が単位行列であることが前提
    return -torch.mean(torch.abs(kurt(Y)))
def ICA(X, loss_func = kurt_loss):
    X = Variable(torch.from_numpy(X), requires_grad =False).float()
    W = Variable((torch.randn((2, 2))), requires_grad =True)
    W.data = W.data / torch.norm(W.data)
    hist_loss = []
    
    for i in range(1000):
        Y = torch.mm(X, W)
        loss = loss_func(Y)
        hist_loss.append(loss.data.numpy())
        loss.backward()

        W.data = W.data - W.grad.data * 0.1
        W.grad.data.zero_()
        W.data = W.data / torch.norm(W.data)
    return Y, W, hist_loss
Y, W, hist_loss = ICA(X, loss_func=kurt_loss)
plt.plot(Y.data.numpy())
plt.show()
plt.plot(hist_loss)

f:id:ksknw:20180519231733p:plain

f:id:ksknw:20180519231743p:plain

以上より,いい感じに復元することができた.

Negentropyを用いた誤差関数

上ではガウス分布度合いを測る指標として,尖度を用いた. しかし,尖度は外れ値に対して敏感であることが知られている. そこで、ガウス分布度合いを測る指標として,Negentropyを用いる方法が提案されている.

 J(\textbf{y}) = H(\textbf{y}_\mathrm{gauss}) - H(\textbf{y})

ここで Hエントロピーガウス分布エントロピーが最も高い確率分布であるため,その差を使ってガウス分布度合いを評価することができるということ. しかし,結局,この計算は確率密度関数が必要なので,データ点だけから求めるのは一般には難しい. そこで,以下のように近似的にNegentropyを求める.

 J(\textbf{y}) \approx \sum _{i=1}^p k_i { \mathbb{E} [G_i(y) ] - \mathbb{E} [G_i(\nu) ] }^2 期待値をサンプルから近似することで,この式は計算できる. ただし  \nu \sim \mathcal{N}(0,\mathbb{I}), k_i>0 . また, Gには色々な関数が設定できるが,例えば, G_1(u) = 1/a_1 \log \mathrm{cosh} a_1 uがある.

def negentropy_loss_logcosh(Y, a=1):
    return -torch.mean(torch.log(torch.cosh(a * Y))/a )**2 
Y, W, hist_loss = ICA(X, loss_func=negentropy_loss_logcosh)
plt.plot(Y.data.numpy())
plt.show()
plt.plot(hist_loss)

f:id:ksknw:20180519231759p:plain

f:id:ksknw:20180519231815p:plain

おわりに

独立成分分析の理解のために,自分で実装してみた. 今回は,天下り的に独立とガウス分布らしさの関係について考えたが,論文では同じ最適化問題相互情報量を最小化するという考えからも導出されている. また,今回は雑に最適化問題を解いたが,FastICAでは不動点反復法(?)で最適化するらしい. 制約付きの最適化問題なので,リーマン多様体上の最適化としても解ける気がする.

参考

リーマン多様体上の最適化 (シュティーフェル多様体上の最適化まで)

概要

  • リーマン多様体上での最適化のお勉強をしている.
  • 多様体上での最適化では接空間上での勾配を求めたり,多様体上へレトラクションしたりする.
  • 例として球面上で,自明な最適化とレイリー商(対称行列の固有ベクトルを求める)をやる.
  • レイリー商から複数の固有値を同時に求めようとするとシュティーフェル多様体上での最適化になる.
  • シュティーフェル多様体はわかった気がするけど,グラスマン多様体がわからんのでできてない.
  • 実装は色んなライブラリがあるけど,pytorchを使う.autogradがあればなんでもいい.

はじめに

年末からお正月で多様体上の最適化について勉強している. 今回は自分の理解のために勉強したことをまとめるが,多様体がそもそもよくわかってないレベルなので色々許してください.

やりたいのはユークリッド空間上での制約つき最適化. 制約付きの最適化はラグランジュの未定乗数法とかで解くイメージだが,ニューラルネットとかがくっついているとどうやって解くのかわからない(あるのかもしれないが). 簡単な方法としてペナルティ項として制約条件を表し,制約なし最適化として解く方法が考えられるが,制約条件を満たす保証はなかったり,パラメータ調整が難しかったりする.

ユークリッド空間上での制約付きの最適化は,制約条件を満たす全ての点があるリーマン多様体上にあるとき,リーマン多様体上での制約なし問題として解くことができる. 例えば,以下は後で使うレイリー商(を変形した式)である.この制約付きの最適化は半径1の球面上の制約なし最適化問題として解釈できる.

 \mathrm{min}_{\bf x} \ {\bf x}^\top A {\bf x}\ \ s.t. \ \ ||x||_2=1

ただし{\bf x} \in \mathbb{R}^n ,\ \ A \in \mathbb{R}^{n\times n}

普通の勾配法に加わる操作としては以下の2つである.

  1. ユークリッド空間で求めた勾配 \nabla _x\  f(x)多様体の接空間に射影する( \mathrm{grad}\ f(x)).
  2. 更新後の点(x_t + \xi多様体からはみ出す)を多様体上に戻す(レトラクション  R_x (x_t+\xi)).

図にすると以下のような感じ.

f:id:ksknw:20180103100907p:plain

まだあんまり理解できてないので数学的な議論を避けて,具体例をやっていく. 最適化にもニュートン法とか共役勾配法とか色々あるけど,今回は単純な最急降下法だけを考える.(学習率も固定) pytorchを使ってfloatで計算しているので,数値誤差がそれなりにある.

こちらなどを参考にして勉強した(している). 数学的な話を含めてわかりやすいように思うのでおすすめ.

以下の内容は主に,これらに書かれていることを自分の理解した範疇で書きなおして,理解のためにpython (pytorch) で実装したというもの.

またこれが有名らしい(?)けど,あまり読めてない (英語苦手...)

3次元空間中の球面の例

まずは簡単な例として球面上の最適化を考える. 最適化する関数は以下のようなもの

\mathrm{argmax}_{\bf x \in \mathbb{R}^3} f({\bf x})={\bf x}[2] \ \ s.t.\ || {\bf x} ||_2=1

解は自明で{\bf x}=[0,0,1]^\topだが,これを多様体上の最急降下法で解くことにする.

ペナルティ項としての実装

リーマン多様体上での最適化を行う前に比較手法として,制約をペナルティ項として表すことを考える. つまり以下の関数を制約なし最適化問題として解く.

\mathrm{argmin}_{\bf x \in \mathbb{R}^3} - x[2] + \gamma |(1-||x||_2)|

コードは以下.

求まった解は

([  5.31967416e-05,   3.33191572e-05,   1.02169939e+00])

で0.02ぐらいずれていて,制約を満たせていない.

収束途中の遷移は以下. 青い球が制約条件で,オレンジの点が上に向かって進んでいる.

多様体からはみだしたり戻ったりを繰り返しながらガタガタしつつ進んでいく. 最急降下法で最適解に落とすには学習率を減衰させる必要があるが,今回はやってないので,こんな感じになる.

f:id:ksknw:20180102211837p:plain

これでもこの問題に対してはある程度の解は得られるが,今回はリーマン多様体上によって制約条件を常に(数値的には)厳密に満たしつつ最急降下法を行う.

多様体上の最適化

以下の最適化問題をリーマン多様体上での最適化として解く. つまり球面をS^{n-1}とするとき,目的関数を以下のように変形する.

\mathrm{argmax}_{\bf x \in S^{n-1}} f({\bf x})={\bf x}[2]

そもそもリーマン多様体が何かあまり理解できてないが, 多様体上の任意の点の接空間でリーマン計量\lt \cdot, \cdot> (内積的なもの)が定義できていればいい(たぶん). 内積が定義できれば接ベクトルの長さが定義できるので,測地線の長さがわかるということ(なのかな?)

ということで,まずは球面S^{n-1}とするときの内積を以下のように定義する.

\lt \xi,\eta> =(\xi^\top \eta)\ \ \xi,\eta\in T_xS^{n-1}\ \ x\in S^{n-1}

(これはいわゆる普通の内積だが,後で行列多様体を考える練習として順にやる.)

接空間上での勾配を考える. ユークリッド空間上での勾配は[0,0,1]^\topであるが,これは多様体の接空間上にはない.

リーマン多様体上での点pにおける勾配(\mathrm{grad}\ f)_pは以下を満たすものとして定義される.  \xiは任意の接ベクトル.

\lt (\mathrm{grad}\ f)_p, \xi>_p = \xi^\top \nabla f

なので,多様体上での勾配は以下のようになる.

\mathrm{grad}f\left({\bf x}\right)= \left(\mathbb{I}-{\bf xx}^\top  \right) \nabla f

\xi^\top x=0なのでこの勾配は上の定義を満たす.

次にレトラクションを考える.今回は球面上に点を戻せばいいので,レトラクションは単純にノルムでわるだけにする.

 R_x(\eta)=\frac{x+\eta}{||x+\eta||_2}

以上から,勾配法で解くプログラムは以下のようになる.

得られた解は以下. float32なのでこれぐらいの誤差は許してほしいところ.

array([  1.29403008e-04,  -1.01325495e-04,   1.00000000e+00], dtype=float32)

制約条件は途中の点でも全て同様に満たされている.

f:id:ksknw:20180102223542p:plain

レイリー商の例

レイリー商(を変形したもの)である以下の最適化問題を考える.

 \mathrm{min}_{\bf x} \ {\bf x}^\top A {\bf x}\ \ s.t. \ \ ||x||_2=1

これも上と同様にリーマン多様体(球面)上の最適化問題として解ける. ちなみにこの問題の解は Aの最小固有値に対応する固有ベクトルになる.

このとき勾配,レトラクションは上と同様.

\mathrm{grad}f\left({\bf x}\right)= \left(\mathbb{I}-{\bf xx}^\top  \right) \nabla f

 R_x(\eta)=\frac{x+\eta}{||x+\eta||_2}

プログラムは以下.

図中青がペナルティ法での点の軌跡.オレンジが多様体上で最適化した時の点の軌跡. 紫は解(固有ベクトル)を表している.

オレンジは多様体上を移動しているのに対して,青は制約条件を満たさない点を経由して最適解へと近づいていく. また,(学習率などの問題もあるが)今回の問題では多様体上での最適化のほうが圧倒的に収束が速かった.

f:id:ksknw:20180102230356p:plain

レイリー商(シュティーフェル多様体上での最急降下法)

ここまではベクトル{\bf x}の最適化について考えてきた. 次に行列 X\in \mathbb{R}^{n\times m}の最適化について考える.

解く問題は先ほどと同様に対称行列A固有ベクトルを求めるもので以下の最適化問題である.

 \mathrm{min}_{X\in \mathbb{R}^{n\times m}} \ \mathrm{trace} \ X^\top A X\ \ s.t. \ \ X^\top X=\mathbb{I}_n

ここでA固有値n個あるがここではそのうちm個を求めるときを考えている。つまりXは正方行列に限らない。

これを解きたいがトレースは行列を並び替えても値が変わらないのでこのままでは解けない。 そこで以下のBrockett cost functionを導入する.

 \mathrm{min}_{X\in \mathbb{R}^{n\times m}} \ \mathrm{trace} \ X^\top A X N \ \ s.t. \ \ X^\top X=\mathbb{I}_n

N=\mathrm{diag}(\mu_1,...,\mu_n),\ \ 0\lt\mu_1\lt...\lt \mu_n

Nによって重み付けされるので,固有値の大きさによって解が順序付けられた状態のみが最適解になる.

制約条件 X^\top X=\mathbb{I}_nを満たす行列からなる多様体をシュティーフェル多様体\mathrm{St(m,n)}という。 まずは 接空間上のリーマン計量を以下の標準内積で定義する。 \lt \xi,\eta>_X = \mathrm{trace} (\xi^\top \eta) \ \ \xi,\eta \in T_X\mathrm{St}(m,n)

(あんまりわかってないので中略) このとき,多様体上の勾配は以下(たぶん). \mathrm{grad}\ f(X) = \nabla f(X) - X \mathrm{symm}(X^\top  \nabla f(X) ), \ \ \mathrm{symm}(D) = 0.5 (D+D^\top)

次にレトラクションを考える. シュティーフェル多様体のレトラクションとして

  • 極分解
  • qr分解を使う方法

などがあるが今回はqr分解を使う方法を用いる.

以上よりプログラムは以下.

これにより2つの固有値を同時に求めることができた.

f:id:ksknw:20180102231107p:plain

おわりに

リーマン多様体上での最急降下をやった. ペナルティで実装するよりは良さそう. 接空間を求めてからどうやって接空間上の勾配を求めれば良いのかがまだわかっていない. あとはQR分解とかの計算コストがどれぐらいなのかを考えたほうが良いかもしれない.

今回はpytorchで実装したが,pymanopt (python manifold optimization) というプロジェクトがあり,このライブラリを使うとシュティーフェル多様体上での最適化などが簡単にできるっぽいので良さそう.(Tensorflowなどとも連携できる(?))

そのうち中略をちゃんと埋める. また,Brockett cost functionではなく,グラスマン多様体として最適化する方法もやっていきたい.

参考

JuliaとPythonで擬似コードや数式を動かそうとしたときの比較

この記事はJulia Advent Calendar 2017の17日目の記事です.

普段はpythonばかり書いていて,juliaは最近文法覚えてきたかなレベルなので色々許してください.

コードの全体はここにあります.

github.com

概要

  • この記事では擬似コードや数式を可能な限りそのままプログラムすることを目的とします.
  • unicode文字を使いまくって以下の画像のようなプログラムを作成します.
  • juliaとpythonで実装して,書きやすさと実行速度を比較します.
  • 書きやすさが悪化するので型指定はしません.
  • 結論は以下です.
    • juliaのほうが色んなunicode文字が使えるから,書きやすく可読性が高い.
    • インデックスが1から始まるのがいい.
    • juliaのほうが倍程度速くなることもあるけど,思ったより速くならない (型指定してないから)
    • juliaのeinsumを何も考えずに使うと激遅になる.(僕がやらかしてるかも)

unicode文字は以下の埋め込みではおそらく微妙な感じに表示されますが,ちゃんと設定されたエディタやjupyter notebookでは以下のように表示されます. f:id:ksknw:20171205223939p:plain

Table of Contents

はじめに (式or擬似コードに可能な限り近いプログラム)

数式を見て実装するときや論文に書いてある擬似コードを実装するとき,可能な限りそれらに近いプログラムを書くようにしている. 数式/擬似コードとプログラムを近づけるために,{ \displaystyle \alpha}などunicode文字を多用する. このようなプログラムを書くことで,

という利点がある.

例えばNo-U-Turn Sampler Hoffman+, 2011擬似コード(の一部,論文より引用)は

f:id:ksknw:20171130224045p:plain

これに対して書いたpythonコードは以下.

def BuildTree(θ, r, u, v, j, ε):
    if j == 0:
        θd, rd = Leapfrog(x, θ, r, v * ε)
        if np.log(u) <= (L(*θd) - 0.5 * np.dot(rd, rd)):
            Cd_ = [[θd, rd]]
        else:
            Cd_ = []
        sd = int(np.log(u) < (Δ_max + L(*θd) - 0.5 * np.dot(rd, rd)))
        return θd, rd, θd, rd, Cd_, sd
    else:
        θ_minus, r_minus, θ_plus, r_plus, Cd_, sd = BuildTree(θ, r, u, v, j - 1, ε)
        if v == -1:
            θ_minus, r_minus, _, _, Cdd_, sdd = BuildTree(θ_minus, r_minus, u, v, j - 1, ε)
        else:
            _, _, θ_plus, r_plus, Cdd_, sdd = BuildTree(θ_plus, r_plus, u, v, j - 1, ε)
        sd = sdd * sd * int((np.dot(θ_plus - θ_minus, r_minus) >= 0) and (np.dot(θ_plus - θ_minus, r_plus) >= 0))
        Cd_.extend(Cdd_)

        return θ_minus, r_minus, θ_plus, r_plus, Cd_, sd

pythonではギリシャ文字などの一部のunicode文字を変数に使うことができる. しかし,一部の例えば { \displaystyle
\theta ^+
}{ \displaystyle
\nabla
} などの記号は使うことができないため,微妙に見難い感じになってしまっていた. 探してみるとjuliaで同じことをしている人がいて,こちらのほうがだいぶ良さそうだった.

bicycle1885.hatenablog.com

juliaでは↑であげたような文字に加えて { \displaystyle
\hat a
} みたいな修飾文字や不等号の { \displaystyle
\le
} とかを使える. juliaで使えるunicode文字一覧はこちら

juliaすごい.ということで,以下ではこれまでpythonで実装したコードをjuliaで書きなおし,書きやすさと実行速度を比較する.

NUTSは既にやられているようなので

  • Dynamic Time Warping
  • Stochastic Block Model

について実装する.

juliaは型をちゃんと指定するともっと速くなるが,「そのまま実装する」という目的に反するのでやらない.

ちなみにNUTSのpython実装の話はこっち.

ksknw.hatenablog.com

Dynamic Time Warping

Dynamic Time Warping (DTW) は,2つの連続値系列データのいい感じの距離を求めるアルゴリズム動的計画法で解く. pythonで実装したときの記事はこっち ksknw.hatenablog.com

A global averaging method for dynamic time warping, with applications to clustering (DTW自体の論文ではない) の擬似コードを参考にしてプログラムを書く. 擬似コードは以下.

f:id:ksknw:20171130230716p:plain

f:id:ksknw:20171130230720p:plain

(ただしこの擬似コードは多分間違っているので,m[i,j]の遷移前のインデックスをsecondに入れるように変えた.)

python

δ = lambda a,b: (a - b)**2
first = lambda x: x[0]
second = lambda x: x[1]

def minVal(v1, v2, v3):
    if first(v1) <= min(first(v2), first(v3)):
        return v1, 0
    elif first(v2) <= first(v3):
        return v2, 1
    else:
        return v3, 2 

def dtw(A, B):
    S = len(A)
    T = len(B)

    m = [[0 for j in range(T)] for i in range(S)]
    m[0][0] = (δ(A[0],B[0]), (-1,-1))
    for i in range(1,S):
        m[i][0] = (m[i-1][0][0] + δ(A[i], B[0]), (i-1,0))
    for j in range(1,T):
        m[0][j] = (m[0][j-1][0] + δ(A[0], B[j]), (0,j-1))

    for i in range(1,S):
        for j in range(1,T):
            minimum, index = minVal(m[i-1][j], m[i][j-1], m[i-1][j-1])
            indexes = [(i-1,j), (i,j-1), (i-1,j-1)]
            m[i][j] = (first(minimum)+δ(A[i], B[j]), indexes[index])
    return m

擬似コードや数式ではインデックスを1から始めることが多いが,pythonは0からなので,頭の中でインデックスをずらしながらプログラムを書く必要がある. 他にはv1が添字っぽくない程度で,大体そのまま書けた.

julia

δ(a,b) = (a - b)^2
# first(x) = x[1] firstは元からあるのでいらない
second(x) = x[2]

function minVal(v₁, v₂, v₃)
    if first(v₁) ≤ minimum([first(v₂), first(v₃)])
        return v₁, 1
    elseif first(v₂) ≤ first(v₃)
        return v₂, 2
    else
        return v₃, 3
    end
end

function DTW(A, B)
    S = length(A)
    T = length(B)
    m = Matrix(S, T)
    m[1, 1] = [δ(A[1], B[1]), (0,0)]
    for i in 2:S
        m[i,1] = [m[i-1, 1][1] + δ(A[i], B[1]), [i-1, 1]]
    end
    for j in 2:T
        m[1,j] = [m[1, j-1][1] + δ(A[1], B[j]), [1, j-1]]
    end
    for i in 2:S
        for j in 2:T
            min, index = minVal(m[i-1,j], m[i,j-1], m[i-1,j-1])
            indexes = [[i-1, j], [i, j-1], [i-1, j-1]]
            m[i,j] = first(min) + δ(A[i],B[j]), indexes[index]
        end
    end
    return m
end

endがある分pythonより長い.一方でpythonでは使えない { \displaystyle
v_1}とか { \displaystyle
≤} が使えるので,より忠実な感じになっている.

実際に書いてみるとわかるけど,インデックスが擬似コードと同じなのは結構大事で,全然頭を使わずにそのまま書けた.

実行速度比較

juliaはpythonに比べて速いらしいので実行速度を比較した. シェルのtime を使って実行速度を比較した. コードの全体はここにある.

julia dtw_julia.jl  2.62s user 0.30s system 110% cpu 2.641 total
python dtw_python.py  2.76s user 0.11s system 99% cpu 2.873 total

期待していたよりも全然速くならなかった. 実行時間が短くてコンパイルのオーバヘッドが大きいのかなと思ったから,forで10回実行するようにした結果が以下.

julia dtw_julia.jl  21.97s user 0.66s system 101% cpu 22.355 total
python dtw_python.py  25.81s user 0.78s system 99% cpu 26.591 total

多少速い気がするけど,期待としては数十倍とかだったので,いまいち. 型指定していないとこれぐらいなのかな(?)

Stochastic Block Model

Stochastic Block Model (SBM)は非対称関係データのクラスタリング手法. 隠れ変数の事後分布をサンプリングして近似する(周辺化ギブスサンプラー). 今回の例では特にクラスタ割り当て { \displaystyle
z_1
} の事後分布を求めて,サンプリングする部分をやる.

pythonでの実装したときの記事はこっち. ksknw.hatenablog.com

関係データ学習という本に載っているクラスタzの事後確率に関する数式は以下. { \displaystyle
z_{1,i}}, { \displaystyle
z_{2,j}}をサンプリングする。

他の変数がgivenだとした時の、{ \displaystyle
z_{1,i}}の事後確率は、

ここで、

はガンマ関数で、 はパラメータ

この事後分布を求めてサンプリングする部分をpythonとjuliaで比較する. 全部書くと見づらいので一部だけ,コードの全体は同じくここにある.

python

n_pos = np.einsum("ikjl, ij", np.tensordot(z1, z2, axes=0), X)  # n_pos_kl = n_pos[k][l]
n_neg = np.einsum("ikjl, ij", np.tensordot(z1, z2, axes=0), 1 - X)

m1_hat = lambda i: m1 - z1[i]  # m1_hat_k = m1_hat[k]

n_pos_hat = lambda i: n_pos - np.einsum("kjl, j", np.tensordot(z1, z2, axes=0)[i], X[i])
n_neg_hat = lambda i: n_neg - np.einsum("kjl, j", np.tensordot(z1, z2, axes=0)[i], 1 - X[i])

α_1_hat = lambda i: α1 + m1_hat(i)
a_hat = lambda i: a0 + n_pos_hat(i)
b_hat = lambda i: b0 + n_neg_hat(i)

aᵢhat = a_hat(i)
bᵢhat = b_hat(i)

p_z1ᵢ_left = logΓ(aᵢhat + bᵢhat) - logΓ(aᵢhat) - logΓ(bᵢhat)
p_z1ᵢ_right_upper = logΓ(aᵢhat + np.dot(X[i], z2)) + logΓ(bᵢhat + np.dot((1 - X[i]), z2))
p_z1ᵢ_right_lower = logΓ(aᵢhat + bᵢhat + m2)
p_z1ᵢ = (α_1_hat(i) * exp(p_z1ᵢ_left + p_z1ᵢ_right_upper - p_z1ᵢ_right_lower)).prod(axis=1)
p_z1ᵢ = p_z1ᵢ.real
p_z1ᵢ = p_z1ᵢ / p_z1ᵢ.sum()
new_z1.append(onehot(choice(range(nb_k), p=p_z1ᵢ), nb_k))

数式には{ \displaystyle
\hat{a}
}{ \displaystyle
n^+
} などが頻出するが,pythonではこれらの文字を使うことができない. このため,m1_hatやn_posなどの変数名になってしまっている.

julia

@einsum n⁺[k,l] := X[i,j] * 𝕀z₁[i,k] * 𝕀z₂[j,l]
@einsum n⁻[k,l] := (ones(X)[i,j] - X[i,j]) * 𝕀z₁[i,k] * 𝕀z₂[j,l]

m̂₁ = m(𝕀z₁) - transpose(𝕀z₁[i,:])
@einsum Σ⁺x𝕀z₂[i,l] := X[i,j] * 𝕀z₂[j,l]
@einsum Σ⁻x𝕀z₂[i,l] := (ones(X)[i,j] - X[i,j]) * 𝕀z₂[j,l]
@einsum n̂⁺[k,l] := n⁺[k,l] - 𝕀z₁[i,k] * Σ⁺x𝕀z₂[i,l]
@einsum n̂⁻[k,l] := n⁻[k,l] - 𝕀z₁[i,k] * Σ⁻x𝕀z₂[i,l]

α̂₁ = α₁ + m̂₁
â = a₀ + n̂⁺
b̂ = b₀ + n̂⁻

temp⁺ = zeros(â)
temp⁻ = zeros(â)
temp = zeros(â)
for j in 1:size(temp⁺)[1]
    temp⁺[j,:] = Σ⁺x𝕀z₂[i,:]
    temp⁻[j,:] = Σ⁻x𝕀z₂[i,:]
    temp[j,:] = sum(𝕀z₂,1)
end

@einsum p_z₁[k,l] := exp(logΓ(â + b̂)-logΓ(â)-logΓ(b̂)
    + logΓ(â+temp⁺)+logΓ(b̂+temp⁻)-logΓ(â+b̂+temp))[k,l]
p_z₁ = α̂₁ .* transpose(prod(p_z₁, 2))
p_z₁ /= sum(p_z₁)

𝕀z₁[i,:] = onehot(sample(1:K, Weights(p_z₁[:])), K)

pythonに対してjuliaではα̂₁やn̂⁺などをそのまま変数名に使えて便利.

実行速度比較

python sbm_python.py  8.74s user 0.09s system 354% cpu 2.492 total
julia sbm_julia.jl  たぶん60000sぐらい(終わらない...)

(なにか間違っているかもしれないが,) なんとjuliaのほうが圧倒的に遅くなってしまった. というか同じ量の計算をさせると終わらない... 調べてみるとeinsumがだめっぽいので,以下ではforで書き直す. (とはいえさすがに遅すぎるので何かやらかしている気がする.)

einsumをfor文で書き直す

juliaのマクロ(@のやつ)はコンパイル時に別の式に展開されるらしい. einsumのマクロもそれぞれの@einsumに対して多重のfor文を展開しているらしく,しかも単独でもnumpyのeinsumより遅いこともあるらしい

ということで,普通にfor文による実装に変更して,実行速度を比較し直す. コードは以下.

m̂₁ = m(𝕀z₁) - transpose(𝕀z₁[i,:])
n⁺ = zeros(K, K)
n⁻ = zeros(K, K)

Σ⁺x𝕀z₂ = zeros(K)
Σ⁻x𝕀z₂ = zeros(K)
for l in 1:K
    for j in 1:N₂
        Σ⁺x𝕀z₂[l] += X[i,j] * 𝕀z₂[j,l]
        Σ⁻x𝕀z₂[l] += (1 - X[i,j]) * 𝕀z₂[j,l]
        for k in 1:K
            n⁺[k,l] += 𝕀z₁[i,k] * X[i,j] *  𝕀z₂[j,l]
            n⁻[k,l] += (1 - X[i,j]) * 𝕀z₁[i,k] * 𝕀z₂[j,l]
        end
    end
end

n̂⁺ = zeros(K,K)
n̂⁻ = zeros(K,K)
for k in 1:K
    for l in 1:K
        n̂⁺[k,l] += n⁺[k,l] - 𝕀z₁[i,k] * Σ⁺x𝕀z₂[l]
        n̂⁻[k,l] += n⁻[k,l] - 𝕀z₁[i,k] * Σ⁻x𝕀z₂[l]
    end
end
α̂₁ = α₁ + m̂₁
â = a₀ + n̂⁺
b̂ = b₀ + n̂⁻

p_z₁ = α̂₁
for k in 1:K
    for l in 1:K
        p_z₁[k] *= exp(logΓ(â[k,l] + b̂[k,l])-logΓ(â[k,l])-logΓ(b̂[k,l]) \
            + logΓ(â[k,l]+Σ⁺x𝕀z₂[l])+logΓ(b̂[k,l]+Σ⁻x𝕀z₂[l])-logΓ(â[k,l]+b̂[k,l]+sum(𝕀z₂,1)[l]))
    end
end
p_z₁ /= sum(p_z₁)

𝕀z₁[i,:] = onehot(sample(1:K, Weights(p_z₁[:])), K)

for...endのせいで見難く感じるが,これは普段pythonばかり書いていてfor文に拒絶反応があるからかもしれない. 書いている時の感想としては,行列の向きとか気にしなくていいので,einsumと同じ程度には書きやすい. 実際にfor文を全部とっぱらえば数式と近くなるはず.

実行速度比較

python sbm_python.py  8.74s user 0.09s system 354% cpu 2.492 total
julia sbm_julia_2.jl  3.62s user 0.46s system 114% cpu 3.559 total

無事にpythonの半分以下の時間で処理が終わるようになった. 普段python使っていると無意識的にfor文を避けていたのが, juliaなら何も気にせずに普通にfor文書きまくれば良さそう.

おわりに

DTWとSBMという2つのアルゴリズムについて,pythonとjuliaでの実装を比較した.

書きやすさの改善 実行速度の改善
DTW インデックスが1から.数字の添字 25.81s → 21.97s
SBM (einsum) α̂₁やn̂⁺をそのまま書ける 激遅 (やらかしているかも)
SBM (for) 8.74s → 3.62s

pythonをjuliaで書き直すことで

  • 数式の添字をそのまま書くことができ,書きやすさと可読性が向上した.
  • 素直にfor文を使えば実行速度は改善した.
  • juliaのeinsumは激遅だった.(僕がやらかしてるかも)
  • 型の指定をしていないからか,実行速度の改善は限定的だった.

また,やっていて気づいたこととして

  • jupyter notebook上で実行しつつ書くというやり方だと,コンパイル時間が以外と鬱陶しい
  • a\hat はOKだけど,\hat aはだめとか.\leはOKだけど,\leqqはだめとかunicode関連で微妙につまることがあった.

この記事ではやらなかったこととして以下がある.気が向いたらやろうと思っている.

  • juliaとpythonで型指定してそれぞれどれぐらい速くなるか.
  • pythonのopt-einsumとの比較

はてなだと表示が微妙だけど,jupyter とかでは綺麗に表示されて見やすいよ!(大事なこと)

追記 (2017/12/17)

Qiitaにて,juliaのコードについてのコメントをいただきました.ありがとうございます.

  • DTWについて,こちらに書かれているように(Matrixの要素の型推論とベクトルの代わりにタプルを使用)コードを変更すると,数十倍程度高速化します.
  • SBMについて,グローバル変数になっているパラメータにconstを指定するだけで,実行時間が90%程度になります.

qiita.com

もう少し自分でも高速化のために,あれこれやってみたいと思います.

参考

pytorchで Canonical Correlation Analysis (正準相関分析)の実装

はじめに

  • pytorchの練習も兼ねて,Canonical Correlation Analysis (正準相関分析)をpytorchを使って実装する.
  • 本当は分散共分散行列からなる行列の一般化固有値問題を解くが,今回は勾配法で解を求める.
  • pytorchのプログラムが間違っていないことを確認するためにscikit-learnでもやる.

Canonical Correlation Analysis (正準相関分析)

こちらの資料を参考にプログラムを書く.

主成分分析が多次元の値に対して,分散が大きい方向に射影するアルゴリズムなのに対して,正準相関分析では2つの多次元変数を射影先で相関が大きくなるように射影するアルゴリズムである.

多次元のデータ x\in\mathbb{R}^{T\times K_x}y\in\mathbb{R}^{T\times K_y} 間の正準相関を考える. ここでTはデータ数(系列データのときはデータ長),K_x, K_yx_i,y_iの次元を表す.

x,yの射影ベクトルをそれぞれ a\in\mathbb{R}^{K_x} , b\in\mathbb{R}^{K_y} とする. x,yの平均がそれぞれ0のとき,射影後の分散 S_{xx},S_{yy} と共分散 S_{xy} は, S_{xx} = (a^\top x) (a^\top x)^\top=a^\top V_{xx}a , S_{yy} = (b^\top y) (b^\top y)^\top=b^\top V_{yy}bS_{xy} = (a^\top x) (b^\top y)^\top=a^\top V_{xy}b

よって,射影後の相関\rho\rho = \frac{a^\top V_{xy}b}{\sqrt{a^\top V_{xx}a}\sqrt{b^\top V_{yy}b}}になる. ここで,aとかbを単純にk倍しても相関は変化しない. 最適化のときに解が複数あるのは面倒なので,制約条件a^\top V_{xx}a = b^\top V_{yy}b = 1を加える.

以上から以下を解くことで正準相関が求まる. \max_{a,b} a^\top V_{xy}b \ \mathrm{s.t.}\ a^\top V_{xx}a = b^\top V_{yy}b = 1

これはラグランジュの未定乗数法で解くことができて,一般化固有値問題になるが,今回はこの式のまま勾配法で解く.

生成データでのテスト

以下のように適当に生成したデータを使って,正準相関分析のプログラムをテストする.

以下では x,y,\mathrm{dummy}\in\mathbb{R}^{200 \times 2}の3つの信号を使う. x,yはもともと同じ信号に大きめのノイズが乗ったもの.dummyは適当なノイズだけのもの. x,y間では正準相関の値が大きくなり,x,dummy間では正準相関の値が小さくなってほしい.

import pylab as plt
import numpy as np
import torch 
from torch.autograd import Variable
import logging
import coloredlogs
from tqdm import tqdm

logger = logging.getLogger(__name__)
coloredlogs.install(
    #fmt="%(asctime)s %(levelname)s %(message)s"
    fmt="%(levelname)s %(message)s"
)

%matplotlib inline

logger.setLevel("DEBUG")
if logger.level >= logging.WARNING:
    tqdm = lambda x:x
u1 = np.r_[np.zeros(40), np.ones(30), np.ones(50)*2, np.ones(30), np.zeros(20), np.ones(30) ]
u2 = np.r_[np.ones(25), np.zeros(25), np.ones(25)*2, np.ones(55), np.zeros(35), np.ones(35) ]
u = np.c_[u1, u2]
plt.plot(u)
plt.show()

f:id:ksknw:20171213230955p:plain

sx = np.random.random(len(u)) * 10
sy = np.random.random(len(u)) * 10
sd = np.random.random(len(u)) * 10
x = np.c_[
    u1 + sx + np.random.random(len(u)),
    u2 - sx + np.random.random(len(u)),   
]
plt.plot(x)
plt.show()

f:id:ksknw:20171213231026p:plain

y = np.c_[
    u1 + sy + np.random.random(len(u)),
    u2 - sy + np.random.random(len(u)),   
]
plt.plot(y)
plt.show()

f:id:ksknw:20171213231039p:plain

dummy = np.c_[
    sd + np.random.random(len(u)),
    -sd + np.random.random(len(u)),   
]
plt.plot(dummy)
plt.show()

f:id:ksknw:20171213231050p:plain

以上のように,見た目では全部ぐちゃぐちゃなノイズにしか見えない3つの信号を生成した. 見た目では全部バラバラだが,x,yに関しては"いい感じの方向"に写像できればノイズを除去できる. x,yそれぞれだけでは"いい感じの方向"がわからないのでノイズと信号を区別できないが,両方を使って相関が大きくなる方向をいい感じの方向と定義することで, 相関する信号を抽出できる.

Scikit-Learnによる実装

sciikit-learnには既に実装されているので,以下のように簡単に正準相関と射影後の変数を求めることができる.

from sklearn.cross_decomposition import CCA
cca = CCA(n_components=1)

x,y間の正準相関

t_x, t_y = cca.fit_transform(x,y)
plt.scatter(t_x, t_y)
print(np.corrcoef(t_x.T, t_y.T )[0,1])
plt.show()

正準相関の値

0.906462806701

射影した後のデータの分布

f:id:ksknw:20171213231110p:plain

x,\mathrm{dummy}間の正準相関

t_x, t_d = cca.fit_transform(x,dummy)
plt.scatter(t_x, t_d)
print(np.corrcoef(t_x.T, t_y.T )[0,1])
plt.show()

正準相関の値

0.141422242479

射影した後のデータの分布

f:id:ksknw:20171213231121p:plain

これと同じ結果を得ることを目標にpytorchで実装する.

PyTorchによる実装

pytorchにはtorch.eigという固有値を求める関数が用意されているので,これを頑張って使えば解けそうではある. が,諸々の事情からペナルティ関数法で制約条件を加えた誤差関数を設定し,autogradで微分を求めて勾配法で解く.

元の最適化問題はこれだった. \max_{a,b} a^\top V_{xy}b \ \mathrm{s.t.}\ a^\top V_{xx}a = b^\top V_{yy}b = 1

制約条件をペナルティ関数にして,以下のような誤差関数を最適化する. \mathrm{argmin}_{a,b} \mathcal{L} = -a^\top V_{xy}b + \gamma \left| 1 - a^\top V_{xx}a \right| + \gamma \left| 1 - b^\top V_{yy}b \right|

\gammaは適当に0.5にした.小さい値だと制約があまり満たされず,逆に大きい値とかにするとうまく収束しないことが多い.この辺がよくわからないので,「機械学習のための連続最適化」という本を買ったけどまだ読めていない.

optimizerは何でもいいけど,単純にSGDにした.pytorch等で実装してると,なんとなくでoptimizerを選べるのでいいなと思う.

x -= np.mean(x, axis=0)
y -= np.mean(y, axis=0)
dummy -= np.mean(dummy, axis=0)
N = len(x)

def CCA(x,y,
        nb_steps = 100000,
        lr = 0.001,
        γ = 0.5):
    vxx = np.dot(x.transpose(), x)/N
    vyy = np.dot(y.transpose(), y)/N
    vxy = np.dot(x.transpose(), y)/N
    
    X = Variable(torch.from_numpy(x)).float()
    Y = Variable(torch.from_numpy(y)).float()
    
    a = Variable(torch.randn(2,1), requires_grad=True)
    b = Variable(torch.randn(2,1), requires_grad=True)
    
    Vxx = Variable(torch.from_numpy(vxx)).float()
    Vyy = Variable(torch.from_numpy(vyy)).float()
    Vxy = Variable(torch.from_numpy(vxy)).float()
    
    hist = []
    from torch import optim
    optimizer = optim.SGD([a, b], lr = 0.0001)

    for i in tqdm(range(nb_steps)):
        Loss = - torch.mm(torch.mm(torch.t(a), Vxy), b)  \
             + γ * torch.abs(1 - torch.mm(torch.mm(torch.t(a),Vxx),a))  \
             + γ * torch.abs(1 - torch.mm(torch.mm(torch.t(b),Vyy),b)) 
        hist.append(Loss.data.numpy()[0,0])
        optimizer.zero_grad()
        Loss.backward()
        optimizer.step()
        
    x_t = np.dot(a.data.numpy().transpose(), x.transpose())
    y_t = np.dot(b.data.numpy().transpose(), y.transpose())
    cca_score = torch.mm(torch.mm(torch.t(a), Vxy), b)
    
    logger.info("a Vxx a: %s"%torch.mm(torch.mm(torch.t(b),Vyy),b).data.numpy())
    logger.info("b Vyy b: %s"%torch.mm(torch.mm(torch.t(a),Vxx),a).data.numpy())
    return hist, x_t, y_t, cca_score

x,y間の正準相関

hist, x_t, y_t, cca_score = CCA(x,y)
plt.plot(hist)
plt.show()
100%|██████████| 100000/100000 [00:19<00:00, 5232.01it/s]
INFO a Vxx a: [[ 1.00017715]]
INFO b Vyy b: [[ 1.00021267]]

lossの履歴

f:id:ksknw:20171213231145p:plain

plt.scatter(x_t, y_t)
plt.show()
print(cca_score)

正準相関の値

Variable containing:
 0.9066
[torch.FloatTensor of size 1x1]

射影した後のデータの分布

f:id:ksknw:20171213231156p:plain

hist, x_t, y_t, cca_score = CCA(x,dummy)
plt.plot(hist)
plt.show()
100%|██████████| 100000/100000 [00:19<00:00, 5199.69it/s]
INFO a Vxx a: [[ 1.00273335]]
INFO b Vyy b: [[ 0.99813366]]

lossの履歴

f:id:ksknw:20171213231211p:plain

x,\mathrm{dummy}間の正準相関

plt.scatter(x_t, y_t)
plt.show()
cca_score

正準相関の値

Variable containing:
1.00000e-02 *
  9.6734
[torch.FloatTensor of size 1x1]

射影した後のデータの分布

f:id:ksknw:20171213231220p:plain

無事にpytorchによる実装で,ある程度正しそうな結果を得ることができた.

おわりに

pytorchで勾配法を使ったCCAの実装をやった. 制約付き最適化を勾配法で解く良いやり方がいまいちわかってないので,勉強したい.

参考

EmacsでブラウザのテキストボックスにMarkdownを書く (校正もする)

概要

Atomic ChromeEmacsからブラウザのテキストボタンを操作しつつ,textlintをflycheckからよんで文章を校正する. flymdでMarkdownのリアルタイムプレビューもする.

こんな感じになる. f:id:ksknw:20171104122228p:plain

はじめに

ブラウザのテキストボックスに文字を打ち込むのが面倒. ブラウザで操作するテキストボックスには,例えば,はてなブログGithubのissueなどがある. これらの操作が面倒な原因として

などがある.これらの問題を解決する方法として以下をやる.

  • Atomic ChromeEmacsからchromeのテキストボックスを操作する
  • textlintでEmacsから文法チェックする
  • flymdでEmacsからMarkdownのリアルタイムプレビューする

(Emacs最高)

Atomic ChromeEmacsからchromeのテキストボックスを操作する

chromeにはAtomic chromeというプラグインがあり,これを使うことでAtomからテキストボックスを操作できる. このプラグインEmacsから使用できるパッケージが公開されている

.emacs.d/Cask に以下を追記してインストール.

(depends-on "atomic-chrome")

また,設定ファイルのどこかに以下を書いておくと,chrome内のプラグインのボタンを押すと勝手にEmacsにバッファが開くようになる. デフォルトでは"C-x C-s"をしても保存しないようになっているが,ファイルもほしいことが多いので,これは保存のままにしておく.

(require 'atomic-chrome)
(atomic-chrome-start-server)
(bind-key* "C-x C-s" 'save-buffer) 

Emacsで書けばキーバインドだけでなく,表をOrgTblで書いたりもできるので良い. Emacsで書いたものをはてなのテキストボックスに貼り付けてもいいけど,画像のアップロードとかが面倒なので,少しだけ楽. はてなもそうだけど,sharelatexとかもEmacsから触れるので色々良い.

textlintでEmacsから文法チェックする

textlintMarkdownなどのテキストファイルの校正を行うアプリケーション. 本体にルールを追加していくことで,自分の用途にあった校正ができる.

作者の方が日本人であることもあって(?)日本語用のルールも充実していて, 例えば,ですますが混在しないようにするとか1文中に同じ助詞が複数でないようにするとか色々ある.

ルールを個別にインストールしていくこともできるが,上位の概念としてプリセットというのがあり,例えば日本語技術文書向けルールのプリセットなどがあり,かつその中に含まれるルールを個別にオフしたりもできるので,これを使うのが良さそう.

以下のようにインストールする.

$ npm i -g textlint textlint-rule-preset-japanese textlint-rule-preset-jtf-style textlint-rule-preset-ja-technical-writing textlint-spellcheck-tech-word textlint-rule-prh

~/.textlintrc に

{
  "rules": {
    "preset-japanese": {
      "sentence-length": false
    },
    "spellcheck-tech-word": true,
    "preset-ja-technical-writing": {
      "ja-no-mixed-period": {
        "periodMark": ".",
      },
      "sentence-length": false
    },
    "prh": {
     "rulePaths" :["~/prh.yml"]
    }
  }
}

とか書いておくと色々動く. 普通に使うときは以下のような感じで使える.

$ textlint hoge.md

こんな感じになる. f:id:ksknw:20171104121530p:plain

これをEmacsから使う方法としてflycheckのルールにする方法がある. ルールの設定はEmacsからでも↑のやつを勝手によんでくれるので,以下のような感じにする.

(require 'flycheck)
(flycheck-define-checker textlint
  "A linter for prose."
  :command ("textlint" "--format" "unix" source-inplace)
 :error-patterns
  ((warning line-start (file-name) ":" line ":" column ": "
            (id (one-or-more (not (any " "))))
            (message (one-or-more not-newline)
                     (zero-or-more "\n" (any " ") (one-or-more not-newline)))
            line-end))
  :modes (text-mode markdown-mode))

(add-to-list 'flycheck-checkers 'textlint)
(add-hook 'markdown-mode-hook 'flycheck-mode)

こんな感じで校正してくれる.

f:id:ksknw:20171104112919p:plain

flymdでEmacsからMarkdownのリアルタイムプレビューする

もともとMarkdown書くときは使っているのでFlymdを使っている. (できればEmacs内でなんらかのcssが反映されたプレビューを見たい).

以上全部やると,最初のような感じになる.

既知の問題点

  • 一旦ファイルを保存しないとtextlintが動かない.
  • flymdはfirefoxで動くので,chromefirefoxEmacsを開いていて頭が悪い
  • 特に表記ゆれについて,自分でやりがちなものを登録していくともっといい感じになると思う.

参考

CycleGANで普通の木を満開の桜にする

はじめに

定期的に生成系のタスクで遊びたくなる. 今回はCycle GANを使って、普通の木を満開の桜に変換してみることにした。

Cycle GAN

論文はこれ. 中身についてはたくさん解説記事があるので、そちらを参考。

Cycle GANでは2つのドメインの間の写像を学習する。 普通のGANとは異なり(乱数ではなく)、片方のドメインの画像を入力して、もう1方のドメインの画像を出力するGeneratorがある。 また、普通のGANの誤差関数(discriminatorの判別)だけでなく、1度片方のGeneratorに通した後、もう一方のGeneratorに通して再構成誤差を計算し、誤差関数に加える。

pix2pixなどでは対になる画像を用意しないと学習ができないが、CycleGANではそういうのがいらないという利点がある。

実験

ImageNetから桜の画像3000枚と普通の木の画像2500枚をダウンロードした. 画像をざっと見た感じ,桜は木全体だけでなく花だけアップの写真が多くて少し気になった. また,桜と富士山が写っている画像も結構多い.

著者の実装githubで公開されているので,それにダウンロードしたデータを食わせる. gtx780tiで2日ぐらい学習させた. 古いGPUGPUメモリが全然足りないので,オプションで色々設定した。(そろそろ買いたい)

$ python train.py --dataroot ./datasets/tree2sakura --name tree2sakura --model cycle_gan --no_dropout --pool_size 50 --loadSize 128 --fineSize 128

結果

1エポック後

元画像 生成画像 再構成画像
f:id:ksknw:20170827212704p:plain f:id:ksknw:20170827212731p:plain f:id:ksknw:20170827212752p:plain
f:id:ksknw:20170827212712p:plain f:id:ksknw:20170827212723p:plain f:id:ksknw:20170827212800p:plain

50エポック後

元画像 生成画像 再構成画像
f:id:ksknw:20170827212240p:plain f:id:ksknw:20170827212326p:plain f:id:ksknw:20170827212256p:plain
f:id:ksknw:20170827212233p:plain f:id:ksknw:20170827212315p:plain f:id:ksknw:20170827212248p:plain

100エポック後

元画像 生成画像 再構成画像
f:id:ksknw:20170827202729p:plain f:id:ksknw:20170827202725p:plain f:id:ksknw:20170827212050p:plain
f:id:ksknw:20170827202852p:plain f:id:ksknw:20170827202846p:plain f:id:ksknw:20170827212103p:plain

100エポックのネットワークで普通の木を桜に変換したもの

白と黒を変換しがちなのが少し気になるけど,なかなかうまくできたと思う.

良さそうなもの

 元画像 生成画像
f:id:ksknw:20170828195918p:plain f:id:ksknw:20170828195929p:plain
f:id:ksknw:20170828200024p:plain f:id:ksknw:20170828200034p:plain

よくなさそうなもの

 元画像 生成画像
f:id:ksknw:20170828200231p:plain f:id:ksknw:20170828200238p:plain
f:id:ksknw:20170828200318p:plain f:id:ksknw:20170828200326p:plain

空が白っぽいとその部分を黒で塗りつぶす傾向がある。なぜかはわからない。

おわりに

既存のコードを動かしただけだけど、思ったよりもうまく動いた。 が、空が異常に暗くなるなど、できてないやつも多い。体感では半々ぐらい。

せっかく対応する画像がなくても学習できるので、こちらでやられている犬猫変換のように、それなりに離れたドメイン間でやると面白そう。

生成系の研究はレッドオーシャンな感じで自分でやるのは大変そうだけど, 公開されているもので遊ぶのは楽しいので、ちょいちょいやっていきたい.

参考

pythonでNUTSの実装

はじめに

前回 ハミルトニアンモンテカルロ法の実装をやった.

今回は No U-Turn Sampler (NUTS)の実装をやる. 論文を参考にした.

コードはここにもある

github.com

NUTS

ハミルトニアンモンテカルロ (HMC)はパラメータの勾配を利用して, 効率的にMCMCサンプリングを行うことができる手法だった.

HMCの問題点は2つ.

  • 更新ステップ数 { \displaystyle L}を適切に指定しなければいけない.
  • 更新の大きさ { \displaystyle \epsilon} を適切に指定しなければいけない.

NUTSは更新ステップ数{ \displaystyle L}を自動的に決定する手法である. 論文内では{ \displaystyle \epsilon}はdual-averaging (Nesterov 2009)を用いて決定するが,今回は決め打ちにする.

更新ステップ数L

ハミルトニアンモンテカルロでは,正規分布によって発生させた運動量を与えて, { \displaystyle L}ステップの間,点を動かす. 予め決められた{ \displaystyle L}ステップの間,点を動かすので,例えば谷にハマった時などガタガタして無駄な計算をしてしまう.

NUTSでは「Uターンしたら点を動かすのをやめる」という規則でこの無駄な計算をなくす. ただし,単純にUターンしたときに中断したら詳細釣り合い条件を満たさなくなるので,少し工夫する.

具体的には,だいたい

  1. 運動量{ \displaystyle r_0}をランダムに決める.
  2. 時間の向き{ \displaystyle v}を{-1, 1}からランダムに選ぶ.
  3. 選んだ時間の向きの方向に{ \displaystyle 2^j}回,点を移動させる(点の軌跡は全て記憶しておく)
  4. どこかでUターンしていないかを確認する
  5. Uターンしていたら6へ,それ以外は{ \displaystyle j+1}して2へ戻る
  6. これまでに計算した軌跡からランダムに1点選んで,サンプリング結果に加えて,1に戻る

という感じ.

実装

以下のようなデータの平均と分散をサンプリングする.

import numpy as np
from numpy import exp
from copy import deepcopy
import pylab as plt
import seaborn as sns
from scipy.stats import norm, gamma
from tqdm import tqdm
from numpy import random
from matplotlib.animation import FuncAnimation

sns.set_style("white")

true_μ = 3
true_σ = 1
nb_data = 1000

x = np.random.normal(true_μ, true_σ, nb_data)

print(x.mean(), x.std())
sns.kdeplot(x)
plt.show()
2.99882646635 1.02180321215

f:id:ksknw:20170806182124p:plain

NUTSのプログラムは以下. 論文のNaive-NUTSの実装をする. Leapfrogで点を移動させるなどはハミルトニアンモンテカルロと同じ.

def log_dh(μ, σ):
    return np.array([-np.sum(x - μ) / σ**2,
                     len(x) / (2 * σ**2) - np.sum((x - μ)**2) / (2 * σ**4)])

def H(θₜ, p):
    return -norm_lpdf(θₜ[0], θₜ[1]) + 0.5 * np.dot(p, p)

def Leapfrog(x, θ, r, ε):
    θ_d = deepcopy(θ)
    r_d = deepcopy(r)
    r_d -= 0.5 * ε * log_dh(θ_d[0], θ_d[1])
    θ_d[0] = θ_d[0] + ε * r_d[0]
    θ_d[1] = θ_d[1] + ε * r_d[1]
    r_d -= 0.5 * ε * log_dh(θ_d[0], θ_d[1])
    return θ_d, r_d
norm_lpdf = lambda μ, σ: np.sum(norm.logpdf(x, μ, σ))
gamma_lpdf = lambda a: np.sum(gamma.logpdf(x, a))

Δ_max = 1000
ε = 0.05
L = norm_lpdf
M = 100

θ0 = np.array([random.randn(), random.gamma(1)])
list_θₘ = [θ0]

{ \displaystyle 2^j} 回の移動は再帰で実装される.

def BuildTree(θ, r, u, v, j, ε):
    if j == 0:
        θd, rd = Leapfrog(x, θ, r, v * ε)
        if np.log(u) <= (L(*θd) - 0.5 * np.dot(rd, rd)):
            Cd_ = [[θd, rd]]
        else:
            Cd_ = []
        sd = int(np.log(u) < (Δ_max + L(*θd) - 0.5 * np.dot(rd, rd)))
        return θd, rd, θd, rd, Cd_, sd
    else:
        θ_minus, r_minus, θ_plus, r_plus, Cd_, sd = BuildTree(θ, r, u, v, j - 1, ε)
        if v == -1:
            θ_minus, r_minus, _, _, Cdd_, sdd = BuildTree(θ_minus, r_minus, u, v, j - 1, ε)
        else:
            _, _, θ_plus, r_plus, Cdd_, sdd = BuildTree(θ_plus, r_plus, u, v, j - 1, ε)
        sd = sdd * sd * int((np.dot(θ_plus - θ_minus, r_minus) >= 0) and (np.dot(θ_plus - θ_minus, r_plus) >= 0))
        Cd_.extend(Cdd_)

        return θ_minus, r_minus, θ_plus, r_plus, Cd_, sd
hist_L = []
hist_C = []
for m in tqdm(range(M)):
    r0 = random.randn(2)
    u = random.uniform(0, exp(L(*list_θₘ[-1]) - 0.5 * np.dot(r0, r0)))

    θ_minus = deepcopy(list_θₘ[-1])
    θ_plus = deepcopy(list_θₘ[-1])
    r_minus = deepcopy(r0)
    r_plus = deepcopy(r0)
    j = 0
    C = [[deepcopy(list_θₘ[-1]), deepcopy(r0)]]
    s = 1

    while s == 1:
        v = random.choice([-1, 1])
        if v == -1:
            θ_minus, r_minus, _, _, Cd, sd = BuildTree(θ_minus, r_minus, u, v, j, ε)
        else:
            _, _, θ_plus, r_plus, Cd, sd = BuildTree(θ_plus, r_plus, u, v, j, ε)

        if sd == 1:
            C.extend(Cd)
        s = sd * int((np.dot(θ_plus - θ_minus, r_minus) >= 0) and (np.dot(θ_plus - θ_minus, r_plus) >= 0))
        j += 1

    index = random.choice(list(range(len(C))))
    list_θₘ.append(C[index][0])

    hist_L.append(L(C[index][0][0], C[index][0][1]))
    hist_C.append(C)
100%|██████████| 100/100 [00:00<00:00, 760.56it/s]
def plot(list_θₘ, hist_C):
    fig = plt.figure()
    ax = fig.add_subplot(111)

    def update(i):
        fig.canvas.draw()
        ax.cla()
        j = i // 3
        if (i % 3) == 0:
            ax.scatter(np.array(hist_C[j])[:, 0, 0], np.array(hist_C[j])[:, 0, 1], linewidth=0, marker=".")
            ax.plot(list_θₘ[:(j), 0], list_θₘ[:(j), 1], c="gray", linewidth=0.3, alpha=0.4)
            ax.scatter(list_θₘ[:(j), 0], list_θₘ[:(j), 1], c="g", linewidth=0, marker=".")

        elif (i % 3) == 1:
            ax.scatter(np.array(hist_C[j])[:, 0, 0], np.array(hist_C[j])[:, 0, 1], linewidth=0, marker=".")
            ax.plot(list_θₘ[:(j + 1), 0], list_θₘ[:(j + 1), 1], c="gray", linewidth=0.3, alpha=0.4)
            ax.scatter(list_θₘ[:(j + 1), 0], list_θₘ[:(j + 1), 1], c="g", linewidth=0, marker=".")
        else:
            ax.scatter(np.array(hist_C[j])[:, 0, 0], np.array(hist_C[j])[:, 0, 1], linewidth=0, marker=".", c="w")
            ax.plot(list_θₘ[:(j + 1), 0], list_θₘ[:(j + 1), 1], c="gray", linewidth=0.3, alpha=0.4)
            ax.scatter(list_θₘ[:(j + 1), 0], list_θₘ[:(j + 1), 1], c="g", linewidth=0, marker=".")

        plot_lim = 30
        if(j) > plot_lim:
            temp_list = np.array([[np.min(np.array(C)[:, 0], axis=0),
                                   np.max(np.array(C)[:, 0], axis=0)] for C in hist_C[j - plot_lim: j]])

            temp_xlim = [np.min(temp_list[:, :, 0]), np.max(temp_list[:, :, 0])]
            xlim_range = temp_xlim[1] - temp_xlim[0]
            temp_ylim = [np.min(temp_list[:, :, 1]), np.max(temp_list[:, :, 1])]
            ylim_range = temp_ylim[1] - temp_ylim[0]
            plt.xlim(temp_xlim[0] - xlim_range * 0.1, temp_xlim[1] + xlim_range * 0.1)
            plt.ylim(temp_ylim[0] - ylim_range * 0.1, temp_ylim[1] + ylim_range * 0.1)

    ani = FuncAnimation(fig, update, frames=len(hist_C) * 3 - 2)
    ani.save("temp.gif", writer="imagemagick", fps=3)
    
plot(list_θₘ, hist_C)

たぶんできていると思う きれいに書き直すつもりだったけど、面倒だったからやめた

「2017年3月14日をもって動画のアップロード機能は終了いたしました」 じゃないんだよなぁ

参考