pythonでガウス過程の実装

はじめに

ガウス過程と機械学習を読んだ. この記事では自分の理解のために,ガウス過程を実装する. notebookはこっちにある.

github.com

今のところの理解では,$N_1$個の入力と出力のペア$(x_1^{(1)}, y_1^{(1)}),\dots ,(x_{1}^{(N_1)}, y_{1}^{(N_1)})$と,出力を求めたい$N_2$個の入力$(x_2^{(1)}, y_2^{(1)}),\dots ,(x_{2}^{(N_2)}, y_{2}^{(N_2)})$があるとき,ガウス過程では以下のようにして出力$y_2$の分布を求める.

  1. 観測$y_1\in \mathbb R^{N_1}$および$y_2 \in \mathbb R^{N_2}$がある1つの多変量ガウス分布に従うとする.$y = [y_1, y_2]^\top$, $x = [x_1, x_2]^\top$とするとき,$y(\mu(x), S(x))$ $y = \mathcal{N}(\mu(x), S(x))$.ここで,平均$\mu$および分散$S$は$x$の関数1
  2. $y_1$が得られたときの,$y_2$の条件付き確率分布を求める.

以下では,平均気温の変動を例に実際にガウス過程を書いてみる.

気温変動の例

この記事では例として以下のような愛知県の気温の変動について考える. データは気象庁のHPから入手できる.

はじめに簡単な例として,ある日の気温$y_1$が観測されたときに,次の日の気温$y_2$の確率分布を求める問題を考える. とりあえず,$y_1$と$y_2$の関係をプロットすると以下.

%matplotlib inline
import numpy as np
import pylab as plt
import pandas as pd
import seaborn as sns
from scipy.stats import multivariate_normal
from random import shuffle

from tqdm import tqdm_notebook as tqdm
data = pd.read_csv("./data.csv")
data.head()
年月日 平均気温 Unnamed: 2 Unnamed: 3
0 2016/8/1 27.9 8 1
1 2016/8/2 26.5 8 1
2 2016/8/3 27.9 8 1
3 2016/8/4 29.0 8 1
4 2016/8/5 29.3 8 1
data = pd.read_csv("./data.csv")
data["年月日"] = pd.to_datetime(data["年月日"])
data = data[data["年月日"] < "2018/8/2"]
temperatures = data.values[:, 1].astype(np.float32)
plt.figure(figsize=(6,6))

plt.scatter(temperatures[:-1],temperatures[1:],  alpha=0.3, linewidths=0, marker=".", c="gray")
plt.xlabel("$y_1$")
plt.ylabel("$Y_2$")
plt.show()

f:id:ksknw:20190813190023p:plain

散布図よりある日の気温と次の日の気温が相関していることがわかる. また,全体的に気温の分布はガウス分布っぽくなっている(ということにする). 以上から,ある日と次の日の気温の分布は多変量ガウス分布$p(y_1, y_2) = \mathcal{N}(\bar{y}, S)$によってモデル化することができそうである. ここで,$\bar{y} \in \mathbb R^{2}$は年間の気温の平均$\bar{y}$を並べたベクトルであり,$S$は分散共分散行列である.

mean_temperature = np.mean(temperatures)
cov_temperatures = np.cov(np.c_[temperatures[:-1], temperatures[1:]], rowvar=0)
Sigma11, Sigma12, Sigma21, Sigma22 = cov_temperatures.reshape(-1)
display(mean_temperature)
display(cov_temperatures)
16.492613

array([[76.93262956, 75.43887714],
       [75.43887714, 77.04298053]])

データより平均は16ぐらい. 分散共分散行列の成分$S_{11}, S_{22}$はともに年間の気温の分散であり,データより$S_{11}, S_{22} = 77$ぐらいで,共分散$S_{12} = S_{21} = 75$ぐらいであった.

データを多変量ガウス分布でモデル化するのが妥当かをなんとなく見るために,この多変量ガウス分布とデータの分布を比較すると以下のようになる.

x = np.linspace(np.min(temperatures)-5,np.max(temperatures)+5)
y = np.linspace(np.min(temperatures)-5,np.max(temperatures)+5)
XX, YY = np.meshgrid(x,y)

shape = XX.shape
XX = XX.reshape(-1)
YY = YY.reshape(-1)
plt.figure(figsize=(12,6))

plt.subplot(121)
plt.title("data distribution")
sns.kdeplot(temperatures[1:], temperatures[:-1], shade=True, cmap="viridis")
plt.scatter(temperatures[1:], temperatures[:-1], 
            alpha=0.3, linewidths=0, marker=".", c="gray")
plt.subplot(122)
plt.title("multivariate gaussian")
plt.contourf(XX.reshape(shape), YY.reshape(shape), 
            multivariate_normal.pdf(np.array(list(zip(XX, YY))), 
                                    mean=[mean_temperature,mean_temperature], 
                                    cov=cov_temperatures).reshape(shape))
plt.scatter(temperatures[1:], temperatures[:-1], 
            alpha=0.3, linewidths=0, marker=".", c="gray")
plt.show()

f:id:ksknw:20190813190102p:plain

(日本の気温は夏と冬で二極化しているので,ガウス分布ではない気もするが) 今回はこれで良しとする.

多変量ガウス分布の条件付き確率

ある日の気温$y_1$と次の日の気温$y_2$の同時分布をモデル化することはできた. 次に,ある日の気温$y_1$が観測されたときの次の日の気温$y_2$の確率,つまり,条件付き確率$p(y_1| y_2) $を考える.

$y_1$が観測されたとき,$y_2$の条件付き確率は $p(y_2 | y_1) = \mathcal{N}(S_{21}S_{11}^{-1}(y_1-\bar{y}) + \bar{y}, S_{22} - S_{21}S_{11^{-1}}S_{12})$ で与えられる.

from scipy.stats import norm
def conditional_pdf(y1, mean, range_y2=np.arange(-3, 3, 0.1)):
    mean = Sigma21/Sigma11*y1 + mean
    variance = Sigma22 - Sigma21 / Sigma11 * Sigma12
    
    return range_y2, norm.pdf(range_y2, loc=mean, scale=variance) * 30

プロットでこれが結局どのような分布なのかを確認する. 例として,今日の気温が25℃,15℃,5℃であったときの$p(y_2|y_1)$を以下に示す.

plt.figure(figsize=(6,6))

plt.scatter(temperatures[1:], temperatures[:-1], 
            alpha=0.3, linewidths=0, marker=".", c="gray")
plt.contour(XX.reshape(shape), YY.reshape(shape), 
            multivariate_normal.pdf(np.array(list(zip(XX, YY))), 
                                    mean=[mean_temperature,mean_temperature], 
                                    cov=cov_temperatures).reshape(shape))

for i in np.linspace(5, 25, 3):
    range_y2, pdf = conditional_pdf(i - mean_temperature, 
                                    mean_temperature,
                                    range_y2=np.linspace(np.min(temperatures)-5, np.max(temperatures)+5,100))
    plt.plot(pdf+i, range_y2, c="C0", alpha=0.8)
plt.axis("equal")
plt.xlabel("$y_1$")
plt.ylabel("$y_2$")
plt.show()

f:id:ksknw:20190813190148p:plain

ガウス分布では$S_{11}$と$S_{21}$が近い値になっているため,$y_1$がある値をとったとき,$y_2$はその値の近くを中心とするガウス分布になることがわかる.

こんな感じで,観測できた変数$y_1$と知りたい値$y_2$(今回は重複していているけど)を1つの多変量ガウス分布でモデル化することができれば,条件付き確率を求めることで,知りたい値$y_2$の分布を求めることができる.

ガウス過程

ここまではある日の気温から次の日の気温を予測する問題を考えてきた. 次にこの問題を少し一般化して,観測できた複数の日の気温$y_1$を使って,期間中の全ての日の気温$y_2$をモデル化する問題を考える.

ids = data.index.values
shuffle(ids)
observed_ids = np.array(sorted(ids[:400]))

data = data.iloc[observed_ids].sort_values("年月日")
y1 = data.values[:, 1].astype(np.float32)
mean_temp = y1.mean()
x1 = observed_ids
x2 = list(set(ids) - set(x1))
xs = np.r_[x1, x2]
data.head()
年月日 平均気温 Unnamed: 2 Unnamed: 3
186 2016-08-02 26.5 8 1
517 2016-08-03 27.9 8 1
144 2016-08-04 29.0 8 1
527 2016-08-05 29.3 8 1
422 2016-08-06 29.8 8 1

観測,予測したい日が複数であった場合にも先ほどと同様に1つの多変量ガウス分布でこれらの同時分布をモデル化する. 今回は$y_2$の平均は$y_1$の平均と同じということにする.

先程の例ではデータから共分散の値を求めていたが,今回は以下のような行列を使ってモデル化する. $$S(x) = K(x; \theta_1, \theta_2) + \theta_3 \mathbb I$$ ここで,$x$は入力,今回の場合は日付(2016/8/1を0とした)である. $K(x; \theta_1, \theta_2)$はカーネル行列であり$K_{nn'}(x_n, x_{n'}; \theta_1, \theta_2) = k(x_n, x_{n'};\theta_1, \theta_2)$,また,今回はガウスカーネル$k(x_n, x_{n'};\theta_1, \theta_2) = \theta_1 \mathrm{exp}\left(- \frac{|x_n - x_{n'}|}{\theta2} \right)$を用いる. $\theta3$は観測のノイズを表す2

import torch
from torch import inverse
y1 = torch.from_numpy(y1)
def gauss_kernel(x, xd, theta1, theta2):
    k = theta1 * torch.exp(- (x-xd)**2/ theta2)    
    return k

def get_covariance(xs, theta1, theta2, theta3):
    cov = []
    for i in xs:
        k = gauss_kernel(
                i,
                torch.Tensor(list(xs)).float(), 
                theta1, theta2
        )
        cov.append(k)
    cov = torch.stack(cov) 
    cov = cov + torch.eye(cov.shape[0]) * theta3
    return cov

適当なパラメータで作った共分散行列は以下. 左上の部分($S_{11}$)は$y_1$に関する分散共分散を表している. 日付順に並んでいるので対角線,つまり,近い日付との共分散が大きくなっていることがわかる. また,右下の部分($S_{22}$)は$y_2$に関する共分散を表しており,ここも同じ傾向がある. 右上($S_{12}$),左下の部分($S_{21}$)は$y_1$,$y_2$の間の共分散を表している. ランダムに間引いて$y_1$と$y_2$を分けているので,日付が近い部分がこういう感じで現れている.

cov = get_covariance(xs=xs, theta1=1, theta2=1, theta3=1)
sns.heatmap(cov, square=True, cmap="viridis", cbar=None)
plt.vlines(x1.shape[0]-1, ymin=0, ymax=xs.shape[0], lw=1, colors="gray")
plt.hlines(x1.shape[0]-1, xmin=0, xmax=xs.shape[0], lw=1, colors="gray")

plt.xlim(0, xs.shape[0])
plt.ylim(xs.shape[0], 0)
plt.show()

f:id:ksknw:20190813190346p:plain

この共分散行列を用いて条件付き確率 $p(y_2 | y_1) = \mathcal{N}(S_{21}S_{11}^{-1}(y_1-\bar{y}) + \bar{y}, S_{22} - S_{21}S_{11^{-1}}S_{12})$を求める.

def conditional_dist(theta1, theta2, theta3, x1, y1, x2, y2=None):
    xs = np.r_[x1, x2]
    y1_mean = torch.mean(y1)
      
    cov = get_covariance(xs, theta1, theta2, theta3 )

    y1_indexes = list(range(len(x1)))
    y2_indexes = list(range(len(x1), len(x1) + len(x2)))
        
    Sigma11 = cov[np.meshgrid(y1_indexes, y1_indexes)]
    Sigma22 = cov[np.meshgrid(y2_indexes, y2_indexes)]
    Sigma12 = cov[np.meshgrid(y2_indexes, y1_indexes)]
    Sigma21 = cov[np.meshgrid(y1_indexes, y2_indexes)]     

    mean = y1_mean + Sigma21 @ inverse(Sigma11) @ (y1-y1_mean)
    variance = Sigma22 - Sigma21 @ inverse(Sigma11) @ Sigma12
        
    if y2 is None:
        return mean, variance
    
    norm = torch.distributions.MultivariateNormal(loc=mean, covariance_matrix=variance)
    return mean, variance, norm.log_prob(y2)
def data4plot(x1, x2, y1, theta1, theta2, theta3):
    xs = np.r_[x1, x2]
    means = []
    variances = []
    plot_range_y2 = torch.stack([torch.linspace(torch.min(y1)-5, 
                                            torch.max(y1)+5, 1000)]).t()
    list_log_probs = []
    for i in tqdm(range(len(xs))):
        mean, variance, log_prob = conditional_dist(theta1=theta1, theta2=theta2, 
                                                    theta3=theta3,
                                                x1=np.array(x1), y1=y1,
                                                x2=np.array([i]), y2=plot_range_y2)
        means.append(mean.data.numpy())
        variances.append(variance.data.numpy())
        list_log_probs.append(log_prob)
    means = np.array(means).reshape(-1)
    variances = np.array(variances).reshape(-1)
    return means, variances, torch.stack(list_log_probs).t(), plot_range_y2
means, variances, log_probs, plot_range_y2 = data4plot(x1, x2, y1, theta1=1, theta2=1, theta3=1)
ids = list(range(len(x1) + len(x2)))
XX, YY = np.meshgrid(ids, plot_range_y2.data.numpy()[:,0])

plt.contourf(XX, YY, np.exp(log_probs), alpha=1)
plt.plot(x1, y1.data.numpy(), c="C1", marker=".")
plt.show()

f:id:ksknw:20190813190452p:plain

なんかだめそうである.一部を拡大してみると以下.

plt.contourf(XX, YY, np.exp(log_probs), alpha=1)
plt.plot(x1, y1.data.numpy(), c="C1", marker=".")
plt.xlim(60, 120)
plt.show()

plt.contourf(XX, YY, np.exp(log_probs), alpha=1)
plt.plot(x1, y1.data.numpy(), c="C1", marker=".")
plt.xlim(160, 220)
plt.show()

f:id:ksknw:20190813190516p:plain

f:id:ksknw:20190813190506p:plain

適当に作ったパラメータでは観測データと結構ずれていてだめそうである.

カーネルの学習

上では共分散行列をカーネルの形で与えることで,観測データを1つの多変量ガウス分布として表し,各日$x_i$での平均気温の確率$p(y_2|y_1)$を求めた. しかし,適当に与えたカーネルのパラメータでは平均気温の変動を正しく捉えることができていない部分があった. これはカーネルのパラメータ$\theta_1$, $\theta_2$,$\theta_3$を適切な値に設定できていないことが原因であると考えられる. そこで,これらを観測データの対数尤度$\log p(y_1 | x_1;\theta)$を最大化するように勾配法を用いて更新する. 参考にした本ではL-BFGSとかを使うとかいてあるが,今回は簡単にAdamを使ってパラメータを更新した. また,$\theta$は全て0以上の値なので,$\theta_i = \exp(η_i), i=1,2,3$として$η_i$を更新することにした.

from torch.optim import LBFGS, Adam
from random import shuffle

eta1 = torch.Tensor([1]).float()
eta2 = torch.Tensor([1]).float()
eta3 = torch.Tensor([1]).float()

eta1.requires_grad = True
eta2.requires_grad = True
eta3.requires_grad = True

optimizer = Adam([eta1, eta2, eta3], lr=0.01)
hist_eta1 = []
hist_eta2 = []
hist_eta3 = []
hist_log_probs = []

for i in tqdm(range(1000)):
    if eta3 < -5:
        eta3.data = torch.Tensor([-5])
    
    theta1 = torch.exp(eta1)
    theta2 = torch.exp(eta2)
    theta3 = torch.exp(eta3)
    
    hist_eta1.append(float(eta1.data.numpy()))
    hist_eta2.append(float(eta2.data.numpy()))
    hist_eta3.append(float(eta3.data.numpy()))
        
    
    means, variances, log_probs = conditional_dist(theta1, theta2, theta3,
                                    x1, y1, 
                                    x1, y1)

    optimizer.zero_grad()
    (-log_probs).backward()
    optimizer.step()
    
    if i==900:
        for param_group in optimizer.param_groups:
            param_group['lr'] *= 0.1

    hist_log_probs.append(log_probs.data.numpy())
plt.plot(hist_eta1)
plt.show()
plt.plot(hist_eta2)
plt.show()
plt.plot(hist_eta3)
plt.show()

f:id:ksknw:20190813190620p:plain

f:id:ksknw:20190813190632p:plain

f:id:ksknw:20190813190642p:plain

事前分布もおかずに最尤推定しているので,$\theta_3$ (観測のノイズの係数)が0に落ちていっている. 最尤推定で求めたパラメータを用いて,$y_2$の平均を求めると以下.

means, variances, log_probs, plot_range_y2 = data4plot(x1, x2, y1, 
                                                       theta1=theta1, theta2=theta2, 
                                                       theta3=theta3)
plt.plot(means)
plt.plot(x1, y1, marker=".")
plt.show()

f:id:ksknw:20190813190704p:plain

先ほどのパラメータよりもできてそうである.

また,$y_2$の分布をプロットすると以下(観測がある点と観測がない点で値が違いすぎるのでheatmapでは書きにくかった). 薄い色が2σ,濃い色が1σの領域を表す. 一部を拡大したものをその下に示す.

plt.plot(means)
plt.fill_between(range(len(means)), means-variances**0.5, means+variances**0.5, facecolor="C0", alpha=0.2)
plt.fill_between(range(len(means)), means-variances**0.5*2, means+variances**0.5*2, facecolor="C0", alpha=0.2)
plt.plot(x1, y1, marker=".")
plt.show()

plt.plot(means)
plt.fill_between(range(len(means)), means-variances**0.5, means+variances**0.5, facecolor="C0", alpha=0.2)
plt.fill_between(range(len(means)), means-variances**0.5*2, means+variances**0.5*2, facecolor="C0", alpha=0.2)
plt.plot(x1, y1, marker=".")

plt.xlim(60, 120)
plt.show()

plt.plot(means)
plt.fill_between(range(len(means)), means-variances**0.5, means+variances**0.5, facecolor="C0", alpha=0.2)
plt.fill_between(range(len(means)), means-variances**0.5*2, means+variances**0.5*2, facecolor="C0", alpha=0.2)
plt.plot(x1, y1, marker=".")

plt.xlim(160, 220)
plt.show()

f:id:ksknw:20190813190724p:plain

f:id:ksknw:20190813190733p:plain

f:id:ksknw:20190813190744p:plain

観測が存在した点では観測点の一点のみ高い確率になり,観測が存在しなかった点では分布が広がっていることがわかる. 今回のように最尤推定すると,$\theta_3$を0に近づけていくので,今回のような分布に収束する. 試しに$\theta_3=0.1$で固定してみると以下のようになる.

from torch.optim import LBFGS, Adam
from random import shuffle

eta1 = torch.Tensor([1]).float()
eta2 = torch.Tensor([1]).float()
eta3 = torch.Tensor([-1]).float()

eta1.requires_grad = True
eta2.requires_grad = True

optimizer = Adam([eta1, eta2], lr=0.01)

hist_eta1 = []
hist_eta2 = []
hist_eta3 = []
hist_log_probs = []

for i in tqdm(range(2000)):
    theta1 = torch.exp(eta1)
    theta2 = torch.exp(eta2)
    theta3 = torch.exp(eta3)
    
    hist_eta1.append(float(eta1.data.numpy()))
    hist_eta2.append(float(eta2.data.numpy()))
    hist_eta3.append(float(eta3.data.numpy()))
        
    
    means, variances, log_probs = conditional_dist(theta1, theta2, theta3,
                                    x1, y1, 
                                    x1, y1)

    optimizer.zero_grad()
    (-log_probs).backward()
    optimizer.step()
    
    if i==900:
        for param_group in optimizer.param_groups:
            param_group['lr'] *= 0.1

    hist_log_probs.append(log_probs.data.numpy())
plt.plot(hist_eta1)
plt.show()
plt.plot(hist_eta2)
plt.show()
plt.plot(hist_eta3)
plt.show()

f:id:ksknw:20190813190942p:plain

f:id:ksknw:20190813190950p:plain

f:id:ksknw:20190813190959p:plain

means, variances, log_probs, plot_range_y2 = data4plot(x1, x2, y1, 
                                                       theta1=theta1, theta2=theta2, 
                                                       theta3=theta3)
plt.plot(means)
plt.plot(x1, y1, marker=".")
plt.show()

f:id:ksknw:20190813191024p:plain

plt.plot(means)
plt.fill_between(range(len(means)), means-variances**0.5, means+variances**0.5, facecolor="C0", alpha=0.2)
plt.fill_between(range(len(means)), means-variances**0.5*2, means+variances**0.5*2, facecolor="C0", alpha=0.2)
plt.plot(x1, y1, marker=".")
plt.show()

plt.plot(means)
plt.fill_between(range(len(means)), means-variances**0.5, means+variances**0.5, facecolor="C0", alpha=0.2)
plt.fill_between(range(len(means)), means-variances**0.5*2, means+variances**0.5*2, facecolor="C0", alpha=0.2)
plt.plot(x1, y1, marker=".")

plt.xlim(60, 120)
plt.show()

plt.plot(means)
plt.fill_between(range(len(means)), means-variances**0.5, means+variances**0.5, facecolor="C0", alpha=0.2)
plt.fill_between(range(len(means)), means-variances**0.5*2, means+variances**0.5*2, facecolor="C0", alpha=0.2)
plt.plot(x1, y1, marker=".")

plt.xlim(160, 220)
plt.show()

f:id:ksknw:20190813191036p:plain

f:id:ksknw:20190813191055p:plain

f:id:ksknw:20190813191104p:plain

本当はパラメータに事前分布とかおいて事後分布を推定したほうがいいように思えるが,今回はここまで.

おわりに

ガウス過程と機械学習の3章までの内容についてコードを書いた. 無限次元とかよくわからないというか無限が何もわからないのでそのへんが理解できているのか怪しい. GPLVMの実装もそのうちやりたい.

参考

ガウス過程と機械学習


  1. はてなの数式でどうしてもΣが打てなかったので諦めた.

  2. (この項がないと$x_1=x_2$のとき分散が0になって死ぬ).

Robust PCA (Principlal Component Pursuit) の実装

はじめに

多変量データは実世界の様々なところで現れる(e.g., 画像、音声、動画).これらの多変量データの多くは、データ自体がもつ次元 (e.g., ピクセル数)よりも小さい次元(自然画像の多様体的なやつ)で表現することができる. データから低次元の構造(低ランク行列)を取り出す代表的な方法として、PCAがある. しかし、実世界で観測されるデータには外れ値が含まれていることが多く、PCAはこの外れ値の影響を強く受けて低ランク行列を求めてしまうという性質がある. 外れ値にロバストに低ランク行列を推定する方法 (a.k.a., robust PCA)として、Principal Component Pursuitが提案されている. ここではPCPの簡単な説明と実装を行い、PCAと比較する.

低ランク行列

なんらかの低ランクな構造を持っている(日本語あやしい)多変量データ X\in \mathbb{R}^{d \times N}は低ランクな行列 L\in \mathbb{R}^{d \times N}とノイズ E\in \mathbb{R}^{d \times N}の和としてかける.ここで、 Nはデータ数、 dは各データの次元を表す.

 X = L + E

ここでは例として、MNISTのデータを Xとして、低ランクな行列 Lを求める問題を考える.

データのダウンロードなどは以下. 最近scikit-learnのMNIST downloadはうまくいかないらしいので,こちらを参考にしてダウンロードした.(いつのまにかdeprecatedになっていた)

以降では Xの平均は0であるとする. X \boldsymbol{1} = \boldsymbol{0}

import numpy as np
import pylab as plt
from sklearn.decomposition import PCA
from sklearn import datasets
%matplotlib inline
import urllib.request
import os.path

url_base = 'http://yann.lecun.com/exdb/mnist/'
key_file = {
    'train_img':'train-images-idx3-ubyte.gz',
    'train_label':'train-labels-idx1-ubyte.gz',
    'test_img':'t10k-images-idx3-ubyte.gz',
    'test_label':'t10k-labels-idx1-ubyte.gz'
}
dataset_dir = "."

def _download(filename):
    file_path = dataset_dir + '/' + filename
    if os.path.exists(file_path):
        return print('already exist')
    print('Downloading ' + filename + ' ...')
    urllib.request.urlretrieve(url_base + filename, file_path)
    print('Done')

def download_mnist():
    for v in key_file.values():
       _download(v)

download_mnist()
def load_mnist(filename):
    img_size = 28**2
    file_path = dataset_dir + '/' + filename
    with gzip.open(file_path, 'rb') as f:
        data = np.frombuffer(f.read(), np.uint8, offset=16)
    return data.reshape(-1, img_size)
import gzip

nb_samples = 3000

data = load_mnist(key_file['train_img'])[:nb_samples]
data = data/255
img1 = data[0].reshape(28, 28)
import pylab as plt
plt.imshow(img1)

f:id:ksknw:20181117140525p:plain

with gzip.open(key_file["train_label"], 'rb') as f:
    labels = np.frombuffer(f.read(), np.uint8, offset=8)[:nb_samples]
labels.shape
(3000,)
X = data - data.mean(axis=0)
plt.pcolormesh(X)
plt.show()

f:id:ksknw:20181117140750p:plain

PCA

PCAではノイズ Eとして、ガウス分布から発生するノイズを考えて、以下のような最適化問題によって低ランク行列 Lを求める.(実際にPCAすると"分散が最大になる方向"に軸を置き直すが、その部分は今は考えずに、低ランクに近似、つまり次元削減の部分にのみ着目する.)

 \mathrm{argmin}_ L \|X - L \|^2 _F \ \  \mathrm{s.t.} \  \mathrm{rank}(L) \lt p

ここで \| \cdot \|^2 _Fはフロベニウスノルム. つまり、PCAはl2の意味でデータ Xに最も"近い"低ランク行列 Lを求めているといえる. この問題は X特異値分解することで、以下のように Lを求める方法が知られている

 L = U_p \Sigma _p V_p ^\top ここで X = U \Sigma V^\topであり、 \Sigma_pは特異値 \sigma_1, \dots \sigma_dのうち、大きい方から p個を対角要素に並べた行列.  U_p, V_pはそれに対応した行,列を抜き出したもの.

from numpy.linalg import svd
u, Σ , vh = svd(X, full_matrices=False)
L = (u[:,:2] @ np.diag(Σ[:2]) @ vh[:2])
plt.pcolormesh(L)
plt.show()

f:id:ksknw:20181117140806p:plain

ジャギジャギしていたデータがなんか綺麗になっていて,なんとなく二次元ぐらいで表現されているような気がする (適当). ここままだと見にくいので,PCAと同じようにデータを特異ベクトル(特異値が大きいほうから2次元)の方向へ射影して散布図として表現すると以下のようになる.

L_d =  X @ np.transpose(vh)
plt.scatter(L_d[:,0], L_d[:,1], marker=".", c=labels, cmap="jet")
plt.show()

f:id:ksknw:20181117140832p:plain

これは(決まらない符号を除いて,) scikit learnのPCAによって求まるものと全く同じになる. そもそもscikit-learnのPCAはSVDで解かれている

pca = PCA()
pcaed = pca.fit_transform(data)
plt.scatter(pcaed[:,0], pcaed[:,1], marker=".", c=labels, cmap="jet")
plt.show()

f:id:ksknw:20181117140850p:plain

外れ値が含まれるデータに対するPCA

ここでデータに外れ値が含まれている場合を考える. 例として、MNISTに含まれるデータのうち、1ピクセルが他のピクセルよりも大きめのランダムな値であるとする.

# data[0] = np.random.randn(784) * 4
data[0, 0] = 100
X = data - data.mean(axis=0)
u, s, vh = svd(X, full_matrices=False)
L_d_with_noise =  X @ np.transpose(vh)
plt.scatter(L_d_with_noise[:,0], L_d_with_noise[:,1], marker=".", c=labels, cmap="jet")
plt.show()

f:id:ksknw:20181117140912p:plain

このようにPCAでは少ない数の外れ値に大きく影響された射影が得られてしまうことがある.

Principal Component Pursuit

外れ値は現実のデータにおいてよく現れるが、データに含まれる構造を見たいという目的で低ランク近似をしていた場合、少ない数の外れ値に大きく影響された結果は望ましくない.

そこで、外れ値の存在を考慮して低ランク行列 Lを求めるアルゴリズム(Principal Component Pursuit)が提案されている.

ここで外れ値を S\in \mathbb{R}^{d\times N}と書く.外れ値は全体の一部であるとして、 Sはスパースな行列であると仮定する. Principal Component Pursuitでは、行列 Xを以下のように表すことができると考える.

 X = L + S

スパース制約をナイーブに考えると Sのl0ノルムを小さくしつつ、 Lを求める必要があるが、この問題はNP-hardなので難しい. そこで、PCPでは以下の緩和した問題を考える.

 \min |L|_* + \lambda |S|_1, \ \ s.t.\ L+S=X

ここで | \cdot |_*はトレースノルム(特異値の和)であり, | \cdot |_1はl1ノルム. この制約付き最適化問題拡張ラグランジュ法を用いて解くことができる. PCPの論文ではこの緩和した問題を解くことで得られる L, Sがいくつかの仮定のもとで、はじめの最適化問題の解と一致することを示している (証明ちゃんと読んでないのでたぶん). 仮定の1つとして(??)、パラメータ \lambdaが以下の値であるとしている.

 \lambda = \frac{1}{\sqrt{n_{(1)}}}、ここで、  n_{(1)} = \max (d, N)である

以上から、PCPを実装すると以下のようになる.

X = data
X.shape
(3000, 784)
S = 0
Y = 0
μ = 2
λ = 1/(max(X.shape))**0.5
print(λ)
nb_epochs = 100
0.018257418583505537
def shrink(x, τ):
    shape = x.shape
    x = x.reshape(-1)
    temp = np.apply_along_axis(lambda a: np.max(np.c_[np.abs(a)-τ, np.zeros(len(a))], axis=1), 0, x)
    return np.sign(x.reshape(shape)) * temp.reshape(shape)
from tqdm import tqdm

hist_normY = []
for i in tqdm(range(nb_epochs)):
    u, Σ, vh = np.linalg.svd(X - S + Y/μ, full_matrices=False)
    L = np.matmul(u * shrink(Σ, 1/μ), vh)
    S = shrink(X - L + Y/μ, λ/μ)
    Y = Y + μ*(X - L - S)
100%|██████████| 100/100 [00:57<00:00,  1.73it/s]
plt.scatter(S[:,0], S[:,1], marker=".", c=labels)
plt.show()

f:id:ksknw:20181117141112p:plain

L_d_PCP =  (X - S - Y)@ np.transpose(vh)

結果

  • PCPで求めた低ランク行列
plt.scatter(L_d_PCP[:,0], L_d_PCP[:,1], marker=".", c=labels, cmap="jet")
plt.show()

f:id:ksknw:20181117140940p:plain

  • ノイズがない場合のPCAの結果
plt.scatter(L_d[:,0], L_d[:,1], marker=".", c=labels, cmap="jet")
plt.show()

f:id:ksknw:20181117141007p:plain

  • ノイズがある場合のPCAの結果
plt.scatter(L_d_with_noise[:,0], L_d_with_noise[:,1], marker=".", c=labels, cmap="jet")
plt.show()

f:id:ksknw:20181117141019p:plain

おわりに

外れ値にロバストなPCAとして、Principal Component Pursuitの実装を行った. PCAを行列の低ランク近似だと思うと、スパースな外れ値がのったモデルを考えるのは妥当であるように思える. 実装した結果を見ても、きちんと外れ値を推定することができていた. 今回はMNISTに適当なノイズをのせたデータを用いたが、論文には実際の動画データに適用した例もある(長時間映っている背景部分が低ランク行列で表現され、移動している人が外れ値になる).

この論文はむちゃくちゃ引用されており、Robust Principal Component Pursuitなど(まだロバストにしたいのか)色々発展もありそうなので、チェックしてみたい.

参考

1次元ガウス分布の測地線と双対測地線のプロット

はじめに

情報幾何学の新展開」という本を読んでいる. まだ序盤しか読めてない上に,あまり理解できていないが,自分の理解のために,例として1次元ガウス分布を対象として,以下の導出とプロットをやる.

  • 指数型分布族の標準形および双対座標系
  • ポテンシャル関数,双対ポテンシャル関数
  • 測地線,双対測地線

かなり天下り的にやっている部分が多いので,主にプロットしただけという感もある.

import numpy as np
import pylab as plt
from numpy import pi as π
from numpy import log
from numpy import exp
from mpl_toolkits.mplot3d import Axes3D
e = exp(1)
%matplotlib inline

座標系のプロット

指数型分布の標準形は以下のようなもの.  p({\bf x, \theta}) = \exp\left(\sum_i \theta^{(i)} x_i - \psi({\bf \theta})\right) 1次元ガウス分布 p(x, \mu , \sigma) = \frac{1}{\sqrt{2 \pi \sigma^2} } \exp \left( \frac{(x-\mu)^2}{2 \sigma^2} \right) であるので,これをごちゃごちゃいじって,

 \theta^{(1)} = \frac{\mu}{\sigma^2}, \theta^{(2)} = \frac{1}{2\sigma^2}, x_{(1)}=x, x_{(2)}=-x^2とすると,ごちゃごちゃいじって,

 p(x, \theta) = \exp(\sum_i \theta^{(i)} x_i - \psi({\bf \theta})) - \left( \frac{(\theta^{(1)})^2}{4\theta^{(2)}} + \frac{1}{2}\log \pi -\frac{1}{2} \log \theta^{(2)}\right)

よって,  \psi(\theta) = \frac{(\theta^{(1)})^2}{4\theta^{(2)}} + \frac{1}{2}\log \pi -\frac{1}{2} \log \theta^{(2)} となる.

ここで, \psi(\theta)をプロットすると,以下のように,これが凸関数っぽくなっていることがわかる.(実際にも凸関数である)

def p2θ(μ, σ2):
    θ1 = μ/σ2
    θ2 = 1/(2*σ2)
    return θ1, θ2

def θ2p(θ1, θ2):
    μ = θ1/(2*θ2)
    σ2 = 1/(2*θ2)
    return μ, σ2

def ψ(θ1, θ2):
    return (θ1)**2/ (4*θ2) + 1/2 * log(π) + 1/2*log(θ2) 
def plot_lattice(t1, t2, c, alpha=1):
    for t1_i, t2_i in zip(t1, t2):
        plt.plot(t1_i, t2_i, c=c, alpha=alpha)
    for t1_i, t2_i in zip(t1.transpose(), t2.transpose()):
        plt.plot(t1_i, t2_i, c=c, alpha=alpha)
temp_θ1 = np.linspace(-1,1, 30)
temp_θ2 = np.linspace(1, 3, 30)

θ1, θ2 = np.meshgrid(temp_θ1, temp_θ2)
ψ_θ = ψ(θ1, θ2)
μ, σ2 = θ2p(θ1,θ2)
fig = plt.figure()
ax = fig.add_subplot(111, projection="3d")
ax.plot_wireframe(θ1, θ2, ψ_θ)
plt.xlabel("$\\theta^{(1)}$")
plt.ylabel("$\\theta^{(2)}$")
ax.set_zlabel("$\psi(\\theta)$")
plt.show()

f:id:ksknw:20180814095351p:plain

%matplotlib inline
plt.figure(figsize=(8,4))

plt.subplot(121)
plot_lattice(μ, σ2, "C0", alpha=0.3)
cont = plt.contour(μ, σ2, ψ_θ)
cont.clabel(fmt='%1.1f', fontsize=14)
plt.xlabel("$\mu$")
plt.ylabel("$\sigma^2$")

plt.subplot(122)
plot_lattice(θ1, θ2, "C1", alpha=0.3)
cont = plt.contour(θ1, θ2, ψ_θ)
cont.clabel(fmt='%1.1f', fontsize=14)
plt.xlabel("$\\theta^{(1)}$")
plt.ylabel("$\\theta^{(2)}$")
plt.tight_layout()
plt.show()

f:id:ksknw:20180814095418p:plain

ここで \psi(\theta)微分 {\bf \eta} = \nabla \psi (\theta)を考えると,  \eta^{(1)} = \frac{\theta^{(1)}}{2\theta^{(2)}}, \eta^{(2)} = -\frac{(\theta^{(1)})^2}{(\theta^{(2)})^2}-\frac{1}{2\theta^{(2)}}

凸関数の微分は元の座標系に対して1つ求まり,また,同じものはない. このため, \etaを座標系として使うこともできる.これを双対座標という. これをプロットすると以下.

def θ2η(θ1, θ2):
    η1 = θ1/(2*θ2)
    η2 = -θ1**2/(2*θ2)**2 - 1/(2*θ2)
    return η1, η2

def η2θ(η1, η2 ):
    θ1 = -η1/(η2+η1**2)
    θ2 = -1/2 * 1/(η2+η1**2)
    return θ1, θ2
μ, σ2 = θ2p(θ1,θ2)
η1,η2 = θ2η(θ1, θ2)
plt.figure(figsize=(12,4))
plt.subplot(131)
plot_lattice(μ, σ2, "C0")
plt.xlabel("$\mu$")
plt.ylabel("$\sigma^2$")
plt.subplot(132)
plot_lattice(θ1, θ2, "C1")
plt.xlabel("$\\theta^{(1)}$")
plt.ylabel("$\\theta^{(2)}$")
plt.subplot(133)
plot_lattice(η1, η2, "C2")
plt.xlabel("$\eta^{(1)}$")
plt.ylabel("$\eta^{(2)}$")
plt.tight_layout()
plt.show()

f:id:ksknw:20180814095434p:plain

この座標系でも \psi(\theta)をプロットすると以下.

μ, σ2 = θ2p(θ1,θ2)
η1,η2 = θ2η(θ1, θ2)

plt.figure(figsize=(12,4))

plt.subplot(131)
plot_lattice(μ, σ2, "C0", alpha=0.3)
cont = plt.contour(μ, σ2, ψ_θ)
cont.clabel(fmt='%1.1f', fontsize=14)
plt.xlabel("$\mu$")
plt.ylabel("$\sigma^2$")

plt.subplot(132)
plot_lattice(θ1, θ2, "C1", alpha=0.3)
cont = plt.contour(θ1, θ2, ψ_θ)
cont.clabel(fmt='%1.1f', fontsize=14)
plt.xlabel("$\\theta^{(1)}$")
plt.ylabel("$\\theta^{(2)}$")

plt.subplot(133)
plot_lattice(η1, η2, "C2", alpha=0.3)
cont = plt.contour(η1, η2, ψ_θ)
cont.clabel(fmt='%1.1f', fontsize=14)
plt.xlabel("$\eta^{(1)}$")
plt.ylabel("$\eta^{(2)}$")
plt.tight_layout()
plt.show()

f:id:ksknw:20180814095446p:plain

ここで,双対ポテンシャル関数は  \phi(\eta) = \theta \cdot \eta - \psi \left( \theta(\eta)\right )で与えられるので(?),  \phi(\eta) = -\frac{1}{2} \log(2\pi e) -\frac{1}{2} \log{ -(\eta^{(1)})^2 + \eta^{(2)}} これをプロットして以下.

これも双対ポテンシャル関数も凸関数であるので,その微分を座標系にすることもできる. これを求めると,元の \theta座標系になる.

def ϕ(η1, η2):
    return -1/2*log(2*π*e) - 1/2 * log(-(η1**2 + η2))
temp_η1 = np.linspace(-0.4, 0.4, 30)
temp_η2 = np.linspace(-1, -0.2, 30)

η1, η2 = np.meshgrid(temp_η1, temp_η2)


θ1,θ2 = η2θ(η1, η2)
μ, σ2 = θ2p(θ1,θ2)

ϕ_η = ϕ(η1, η2)
fig = plt.figure()
ax = fig.add_subplot(111, projection="3d")
ax.plot_wireframe(η1, η2, ϕ_η)
plt.xlabel("$\eta^{(1)}$")
plt.ylabel("$\eta^{(2)}$")
ax.set_zlabel("$\phi(\eta)$")
plt.show()

f:id:ksknw:20180814095501p:plain

%matplotlib inline

plt.figure(figsize=(12,4))

plt.subplot(131)
plot_lattice(μ, σ2, "C0", alpha=0.3)
cont = plt.contour(μ, σ2, ϕ_η)
cont.clabel(fmt='%1.1f', fontsize=14)
plt.xlabel("$\mu$")
plt.ylabel("$\sigma^2$")

plt.subplot(132)
plot_lattice(θ1, θ2, "C1", alpha=0.3)
cont = plt.contour(θ1, θ2, ϕ_η)
cont.clabel(fmt='%1.1f', fontsize=14)
plt.xlabel("$\\theta^{(1)}$")
plt.ylabel("$\\theta^{(2)}$")

plt.subplot(133)
plot_lattice(η1, η2, "C2", alpha=0.3)
cont = plt.contour(η1, η2, ϕ_η)
cont.clabel(fmt='%1.1f', fontsize=14)
plt.xlabel("$\eta^{(1)}$")
plt.ylabel("$\eta^{(2)}$")
plt.tight_layout()
plt.show()

f:id:ksknw:20180814095517p:plain

測地線と双対測地線

測地線と双対測地線は,それぞれの座標系で直線として与えられる. つまり測地線は  \theta(t) = t \theta_1 + (1-t)\theta_2, 双対測地線は  \eta(t) = t \eta_1 + (1-t)\eta_2, である.

対応する点  \theta \etaについて,それぞれの測地線をプロットすると以下.

θ_1 = np.array([-1, 1])
θ_2 = np.array([1, 3])

list_t = np.arange(0, 1, 0.01)
line_θ = np.transpose([θ_1*t + θ_2*(1-t) for t in list_t])
line_p1 = θ2p(*line_θ)
line_η = θ2η(*line_θ)

temp_θ1 = np.linspace(-1,1, 30)
temp_θ2 = np.linspace(1, 3, 30)

θ1, θ2 = np.meshgrid(temp_θ1, temp_θ2)
μ, σ2 = θ2p(θ1, θ2)
η1,η2 = θ2η(θ1,θ2)


plt.figure(figsize=(12,4))

plt.subplot(131)
plot_lattice(μ, σ2, "C0", alpha=0.3)
plt.plot(line_p1[0], line_p1[1], "C4")
plt.xlabel("$\mu$")
plt.ylabel("$\sigma^2$")

plt.subplot(132)
plot_lattice(θ1, θ2, "C1", alpha=0.3)
plt.plot(line_θ[0], line_θ[1], "C4")
plt.xlabel("$\\theta^{(1)}$")
plt.ylabel("$\\theta^{(2)}$")

plt.subplot(133)
plot_lattice(η1, η2, "C2", alpha=0.3)
plt.plot(line_η[0], line_η[1], "C4")
plt.xlabel("$\eta^{(1)}$")
plt.ylabel("$\eta^{(2)}$")
plt.tight_layout()
plt.show()

f:id:ksknw:20180814095531p:plain

η_1 = np.array([-0.5, -3/4])
η_2 = np.array([ 1/6, -7/36])

list_t = np.arange(0, 1, 0.01)
line_η= np.transpose([η_1*t + η_2*(1-t) for t in list_t])
line_θ= η2θ(*line_η)
line_p2 = θ2p(*line_θ)
    
temp_η1 = np.linspace(-0.5, 0.5, 30)
temp_η2 = np.linspace(-1, -0.2, 30)

η1, η2 = np.meshgrid(temp_η1, temp_η2)
θ1,θ2 = η2θ(η1,η2)
μ, σ2 = θ2p(θ1, θ2)

    
plt.figure(figsize=(12,4))

plt.subplot(131)
plot_lattice(μ, σ2, "C0", alpha=0.3)
plt.plot(line_p2[0], line_p2[1], "C3")
plt.xlabel("$\mu$")
plt.ylabel("$\sigma^2$")

plt.subplot(132)
plot_lattice(θ1, θ2, "C1", alpha=0.3)
plt.plot(line_θ[0], line_θ[1], "C3")
plt.xlabel("$\\theta^{(1)}$")
plt.ylabel("$\\theta^{(2)}$")

plt.subplot(133)
plot_lattice(η1, η2, "C2", alpha=0.3)
plt.plot(line_η[0], line_η[1], "C3")
plt.xlabel("$\eta^{(1)}$")
plt.ylabel("$\eta^{(2)}$")
plt.tight_layout()
plt.show()

f:id:ksknw:20180814095542p:plain

これらを元のパラメータ空間でプロットすると以下のようになる.

temp_μ = np.linspace(-0.5,0.5, 30)
temp_σ2 = np.linspace(0.01, 1.0, 30)

μ, σ2 = np.meshgrid(temp_μ, temp_σ2)

plot_lattice(μ, σ2, "C0", alpha=0.3)
plt.plot(line_p1[0], line_p1[1], "C4")
plt.plot(line_p2[0], line_p2[1], "C3")

plt.text(-0.2, 0.5,"Dual geodesic line", size=14, color="C3")
plt.text(-0.4, 0.2,"Geodesic line", size=14, color="C4")

plt.xlabel("$\mu$")
plt.ylabel("$\sigma^2$")
plt.tight_layout()
plt.show()

f:id:ksknw:20180814095554p:plain

おわりに

とりあえず双対測地線のプロットまでをやった. 理解できてない部分もあるので,何か間違っているかもしれない.

この本は個人的には読みやすく感じるが,数学に詳しい人に聞くと, 先に接続とか平坦とかをやってから,測地線の方程式を求めたりしたほうがいいらしい (微分幾何を専攻したい人生だった). 実際天下り的に感じる部分もあるので,もうちょっと読めたら,ちゃんとそっちの順番で理解していきたいと思う.

参考

情報幾何学の新展開

機械学習とかの実験のパラメータをいい感じにやるコード

機械学習とかやっていると複数のパラメータを,コマンドラインから設定したいことがある. 割といい感じにできるようになった気がするので,備忘録として書いておく.

コマンドラインパーサとしてclickを使う. clickについては,こちらが詳しくていい感じだった.

この記事では以下ができるようにする.

  1. コマンドラインからいくつかのパラメータを設定する.
  2. パラメータを増やしたくなったときに,さくっとできるようにする.
  3. パラメータを名前に含んだディレクトリを作成して,その中に色々入れる.
import click
import datetime
import os
import pickle


def train(batch_size, learning_rate, nb_epochs, **kwargs):
    return 0, 0


def test(test_mode, **kwargs):
    print(test_mode)
    return 0


@click.command()
@click.option("-b", "--batch_size", type=int, default=12, help="mini batch size")
@click.option("-l", "--learning_rate", type=float, default=0.01)
@click.option("-n", "--nb_epochs", type=int, default=10, help="number of epochs")
@click.option("-m", "--test_mode", type=click.Choice(["hoge", "hogehoge"]))
@click.option("-o", "output_prefix", type=str, default="../results")
def main(**kwargs):
    net, loss_hist = train(**kwargs)
    test(**kwargs)

    output_prefix = kwargs.pop("output_prefix")
    output_dir = output_prefix + "/hoge_"
    now = datetime.datetime.now().strftime("%y%m%d%H%M%S")
    output_dir += now
    kwargs["output_dir"] = output_dir

    os.mkdir(output_dir)
    with open(output_dir + "/log.txt", "w") as f:
        for key, param in kwargs.items():
            print("{}, {}".format(key, param), file=f)
    with open(output_dir + "/params.pkl", "wb") as f:
        pickle.dump(kwargs, f)
    with open(output_dir + "/loss_hist.pkl", "wb") as f:
        pickle.dump(loss_hist, f)

    return


if __name__ == '__main__':
    main()

kwargsには辞書型が入ってくるので,あとはお好みで色々できる. 引数の最後に kwargsをdef train(batch_size, learning_rate, nb_epochs, **kwargs):な感じで入れておけば,必要な変数だけとってきて,残りはおいておくことができるので,いい感じに書けた.

独立成分分析の実装

はじめに

統計的因果探索を読んでいる. この中で,独立成分分析が出てくるが,そういえばあんまり理解できていなかったので,実装して理解する. 日本語の本も見たけど,結局 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

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

参考