部分一致DTW(SPRING)の実装

はじめに

以前,2つの時系列データの距離的なものを測るアルゴリズムであるDTWについて,以下のようなものを書きました.

ksknw.hatenablog.com

DTWはいい感じのアルゴリズムですが,時系列データの全体と全体を対応付けるので,短いパターンを長い時系列データの中から見つけるみたいなタスクに直接用いることはできません(やろうと思えばできるけど計算量やばい). なんかもっといいやつないかなと思っていたところに,たまたま読んだ論文で使われていたSPRINGというアルゴリズムを知ったので,ここでは実装しつつちゃんと理解してみることにしました. 以下ではpythonを使ってSPRINGを実装します.このノートブックの全体は以下にあります.

github.com

目的

ここでは,時系列データ$X$からある時系列パターン$Y$と似ている区間を見つける問題について考えます. 例として,以下の図に示すような時系列データ$X$ (青線)と時系列パターン$Y$(オレンジ線)を考えます. ちなみにこれは気象庁配布されている名古屋の気温データです.

import numpy as np
import pylab as plt
import pandas as pd
import matplotlib.gridspec as gridspec
from tslearn.metrics import dtw_path
data = pd.read_csv("./data.csv", header=None)[1].values
X = data[:1000:4]
Y = data[1000::4]

plt.plot(X, label="X")
plt.plot(Y, label="Y")
plt.legend()
plt.show()

f:id:ksknw:20191228102715p:plain

ここでの目的は以下のような部分時系列(緑部分)を時系列データ$X$から検出することです.

from spring import spring
for path, cost in spring(X, Y, 80):
    plt.plot(X, c="gray", alpha=0.5)
    plt.plot(path[:,0], X[path[:,0]], C="C2")

f:id:ksknw:20191228102754p:plain

2つの時系列データの距離(のようなもの)を測る方法の1つとして,Dynamic Time Warping (DTW)があります. DTWについて詳しくは前書いたやつを参照. しかし,DTWは2つの時系列データ全体同士を比較するものなので,これを今回の目的にそのまま用いることはできません. 時系列データの部分一致の問題を扱うアルゴリズムの1つとして,SPRINGがあります. 以下では,DTWについて簡単に説明したあと,部分一致DTW,そして,SPRINGについて説明します.

DTWと動的計画法

DTWでは2つの時系列データ$X = [X^{(1)}, \dots, X^{(T_x)}]$と時系列データ$Y = [Y^{(1)}, \dots, Y^{(T_y)}]$の距離的なもの(以下ではDTW距離)を求めます. DTWでは2つの時系列データに含まれる各フレーム$X_i$, $Y_j$を,フレーム間の距離の総和が最小になるように対応付けを求めます. このときの,フレーム間の距離の総和がDTW距離です.

ここではフレーム間の距離として二乗距離を考えます. 各フレーム間の距離を並べた行列$\Delta$, $\Delta_{i,j}= \|X^{(i)} - Y^{(j)} \|^2_2$は以下のようなものです. (以下では,説明のために,一旦同じぐらいの長さの時系列データ$X,Y$を考えます.)

X = data[280:400:4]
Y = data[1000::4]
D = (np.array(X).reshape(1, -1) - np.array(Y).reshape(-1, 1))**2
plt.imshow(D, cmap="Blues")
plt.xlabel("$X$")
plt.ylabel("$Y$")
plt.show()

f:id:ksknw:20191228102848p:plain

次にDTWではフレーム間の対応付け(アライメント)を求めます. 実際に求まるアライメントを以下の図に示します.

def dist(x, y):
    return (x - y)**2

def get_min(m0, m1, m2, i, j):
    if m0 < m1:
        if m0 < m2:
            return i - 1, j, m0
        else:
            return i - 1, j - 1, m2
    else:
        if m1 < m2:
            return i, j - 1, m1
        else:
            return i - 1, j - 1, m2

def dtw(x, y):
    Tx = len(x)
    Ty = len(y)

    C = np.zeros((Tx, Ty))
    B = np.zeros((Tx, Ty, 2), int)

    C[0, 0] = dist(x[0], y[0])
    for i in range(Tx):
        C[i, 0] = C[i - 1, 0] + dist(x[i], y[0])
        B[i, 0] = [i-1, 0]

    for j in range(1, Ty):
        C[0, j] = C[0, j - 1] + dist(x[0], y[j])
        B[0, j] = [0, j - 1]

    for i in range(1, Tx):
        for j in range(1, Ty):
            pi, pj, m = get_min(C[i - 1, j],
                                C[i, j - 1],
                                C[i - 1, j - 1],
                                i, j)
            C[i, j] = dist(x[i], y[j]) + m
            B[i, j] = [pi, pj]
    cost = C[-1, -1]
    
    path = [[Tx - 1, Ty - 1]]
    i = Tx - 1
    j = Ty - 1

    while ((B[i, j][0] != 0) or (B[i, j][1] != 0)):
        path.append(B[i, j])
        i, j = B[i, j].astype(int)
    path.append([0,0])
    return np.array(path), cost, C
path, dtw_dist, C = dtw(X, Y)
def plot_path(paths, A, B, D):
    plt.figure(figsize=(5,5))
    gs = gridspec.GridSpec(2, 2,
                       width_ratios=[1,5],
                       height_ratios=[5,1]
                       )
    ax1 = plt.subplot(gs[0])
    ax2 = plt.subplot(gs[1])
    ax4 = plt.subplot(gs[3])

    ax2.pcolor(D, cmap=plt.cm.Blues)
    ax2.get_xaxis().set_ticks([])
    ax2.get_yaxis().set_ticks([])
    
    for path in paths:
        ax2.plot(path[:,0]+0.5, path[:,1]+0.5, c="C3")
    
    ax4.plot(A)
    ax4.set_xlabel("$X$")
    ax1.invert_xaxis()
    ax1.plot(B, range(len(B)), c="C1")
    ax1.set_ylabel("$Y$")

    ax2.set_xlim(0, len(A))
    ax2.set_ylim(0, len(B))
    plt.show()
plot_path([np.array(path)], X, Y, D)

f:id:ksknw:20191228102911p:plain

図の中で,青線,オレンジ線はそれぞれ$X$, $Y$を表しています. 右上の図の青色は$\Delta$を表しており,1本の赤線(パス)で表したものがアライメントされたフレームを表します. ここで,アライメントは以下の3つの制約条件を満たすもののなかで,パス上のコストが最小になるものです.

  • 境界条件: 両端が左下と右上にあること
  • 単調性: 左下から始まり,→,↑,➚のいずれかにしか進まないこと
  • 連続性: 繋がっていること

フレームの対応付けを図に書くと以下です. 境界条件は$X,Y$の両端のフレームを対応付けること,単調性は対応づけがクロスしないこと,連続条件は対応付けられないフレームが存在しないことを表しています.

for line in path:
    plt.plot(line, [X[line[0]], Y[line[1]]], linewidth=0.8, c="gray")
plt.plot(X)
plt.plot(Y)
plt.show()

f:id:ksknw:20191228102926p:plain

DTWではこれらの制約条件を満たすパスの集合$\mathcal A$のうち,パス上の距離$\Delta$の総和が最小となるパスを求めるアルゴリズムです. 式で書くと以下です. $$\min_{A \in \mathcal A} \big < A, \Delta \big >$$ ただし,$\big < \cdot, \cdot \big >$は行列の内積を意味します.

この最適化問題動的計画法を使って解くことができます. 動的計画法では,左下のマスからスタートして,各マスに到達するため最小の累積コストを1マスずつ求めます. このとき,→,↑,➚のいずれかにしか移動できないので,$i+1, j+1$マスの最小の累積和$C$は $$C_{i,j} = \min \big(C_{i, j+1}, C_{i+1, j}, C_{i+1, j+1} \big)$$ で求めることができます. ただし,左下のマス($i=1, j=1$)では,$C_{1,1}=\Delta_{1,1}$であり,下一列($i\neq 1, j=1$)のときは, $$C_{i+1, 1} = C_{i, 1} + \Delta_{i+1, 1}$$ です. また,左一列も同様にして $$C_{i, j+1} = C_{1, j} + \Delta_{1, j+1}$$ です.

累積コストを求めるときに使ったマス($\mathrm{argmin} (C_{i, j+1}, C_{i+1, j}, C_{i+1, j+1})$)を記録しておくことで,終端(右上)にたどり着いた時,この記録を逆にたどることで,累積コストを最小にするパスを求めることができます. この部分のアルゴリズムについては,こちらのアニメーションがわかりやすいです. $C$をプロットすると以下です.

plt.pcolor(C, cmap="Blues")
plt.show()

f:id:ksknw:20191228102940p:plain

DTWを用いることで,時系列データ間の類似度を求めることができますが,今回ような長い時系列データから短い時系列パターンと一致する区間を検出するタスクに直接用いることはできません. ちなみに時系列データと時系列パターンを入力として,DTWをやってみると以下のようなパス,アライメントが得られます.

X = data[:1000:4]
Y = data[1000::4]

D = (np.array(X).reshape(1, -1) - np.array(Y).reshape(-1, 1))**2
path, dtw_dist = dtw_path(X, Y)
plot_path([np.array(path)], X, Y, D)

for line in path:
    plt.plot(line, [X[line[0]], Y[line[1]]], linewidth=0.8, c="gray")
plt.plot(X)
plt.plot(Y)
plt.show()

f:id:ksknw:20191228103002p:plain

f:id:ksknw:20191228103014p:plain

むりやり全体と全体を対応付けようとするため,よくわからない感じになっています.

部分一致DTW

次に,時系列データ$X=[X^{(1)}, \dots, X^{(T_x)}]$から,時系列パターン$Y = [Y^{(1)}, \dots, Y^{(T_y)}], (T_x > T_y)$に,最もDTW距離が近い区間を1箇所見つける問題を考えます. ナイーブな方法として,ありうる全ての部分時系列$[X^{(t)},...,X^{(s)}]$とパターン$Y$とのDTW距離を求め,最も小さいDTW距離をもつ区間を検出する方法が考えられます. しかし,この方法は時系列データ$X$の長さが長くなると計算量が大きくなってしまいます.

SPRINGでは,DTWの動的計画法を一部修正することによって,部分一致問題を解くことを提案しています. (以下では文献とは若干異なるアルゴリズムを説明しますが,たぶんやっていることは同じです.これは自分的実装の楽さから変更されています.) 以下ではこのアルゴリズムを部分一致DTWと呼びます.

DTWの境界条件を以下のように変更することで,部分一致問題に対応することができます.

  • 境界条件: 下の一列のどこかのマスと上の一列のどこかのマスを通る

これは$Y$の両端が$X$のいずれかのフレームと対応付けられることを意味します. 一方で,$X$の両端のフレームは必ずしも対応付けられる必要はありません.

アルゴリズムとしては,動的計画法の中で以下の2点を変更します. - 下一列($i\neq 1, j=1$)の累積コストを$C_{i+1, 1} = C_{i, 1} + \Delta_{i+1, 1}$から$C_{i+1, 1} = \Delta_{i+1, 1}$に変更 - 終了地点を右上のマスから,$\mathrm{argmin}_i (C_{i, T_y})$へと変更

これらの変更はそれぞれ,始点と終点の変更に対応します. これらの変更によって得られたパスを以下の図に示します.

X = data[:1000:4]
Y = data[1000::4]
def dist(x, y):
    return (x - y)**2

def get_min(m0, m1, m2, i, j):
    if m0 < m1:
        if m0 < m2:
            return i - 1, j, m0
        else:
            return i - 1, j - 1, m2
    else:
        if m1 < m2:
            return i, j - 1, m1
        else:
            return i - 1, j - 1, m2

def partial_dtw(x, y):
    Tx = len(x)
    Ty = len(y)

    C = np.zeros((Tx, Ty))
    B = np.zeros((Tx, Ty, 2), int)

    C[0, 0] = dist(x[0], y[0])
    for i in range(Tx):
        C[i, 0] = dist(x[i], y[0])
        B[i, 0] = [0, 0]

    for j in range(1, Ty):
        C[0, j] = C[0, j - 1] + dist(x[0], y[j])
        B[0, j] = [0, j - 1]

    for i in range(1, Tx):
        for j in range(1, Ty):
            pi, pj, m = get_min(C[i - 1, j],
                                C[i, j - 1],
                                C[i - 1, j - 1],
                                i, j)
            C[i, j] = dist(x[i], y[j]) + m
            B[i, j] = [pi, pj]
    t_end = np.argmin(C[:,-1])
    cost = C[t_end, -1]
    
    path = [[t_end, Ty - 1]]
    i = t_end
    j = Ty - 1

    while (B[i, j][0] != 0 or B[i, j][1] != 0):
        path.append(B[i, j])
        i, j = B[i, j].astype(int)
        
    return np.array(path), cost
path, cost = partial_dtw(X, Y)
D = (np.array(X).reshape(1, -1) - np.array(Y).reshape(-1, 1))**2
plot_path([np.array(path)], X, Y, D)

f:id:ksknw:20191228103046p:plain

for line in path:
    plt.plot(line, [X[line[0]], Y[line[1]]], linewidth=0.8, c="gray")
plt.plot(X)
plt.plot(Y)
plt.plot(path[:,0], X[path[:,0]], c="C2")
plt.show()

f:id:ksknw:20191228103117p:plain

$X$から$Y$と似た部分時系列が一部のみ抽出されていることがわかります.

複数の区間の検出

上で説明した手法を用いることで,長い時系列データ$X$から時系列パターン$Y$に一致する部分時系列を抜き出すことができました. 一方で,この方法では最もDTW距離が小さい1つの区間しか検出することができませんでした. SPRINGでは,$X$の中からパターン$Y$とDTW距離が$\epsilon$以下の部分時系列を複数検出することを提案しています. SPRINGでは特にオンライン設定を対象としており,ここでも同様に$X$が徐々に観測されていくときを考えます.

このとき,DTW距離が$\epsilon$以下の部分時系列の区間が重複している場合が考えられます. このような場合に,全ての区間を列挙することもできますが,可視化などの応用を考えると,全ての区間を検出することは適切ではなさそうです. SPRINGでは重複する区間を全て列挙するのではなく,重複するもののうち,最もDTW距離が小さいもの1つを検出します.

SPRINGでは,この問題を扱うために,各マスについて区間の開始地点を記録,伝搬することを提案しています. つまり,各マスを通る最小コストのパスの区間の開始地点$S_{ij}$をそれぞれのマスで記録します. これは動的計画法の中で,以下のように埋めることができます. $$S_{i,1}=i$$ $$S_{i,j} = S_{i^*, j^*}, (i^*, j^*) = \mathrm{argmin}(C_{i,j+1}, C_{i+1,j}, C_{i,j})$$

この操作によって,$S_{i,T_y}$には,$(i, T_y)$を通るコスト最小のパスの開始位置が記録された状態になります. この$S$を用いることで,以下の手順でDTW距離がしきい値$\epsilon$以下の区間を検出することができます. さらに,ここがよく出来ているところですが,SPRINGでは区間を検出するとき,今後観測される$X$に検出した区間と重複し,かつ,DTW距離が検出した区間よりも小さい区間が存在しないことを保証することができます.

  • 現在の時刻を$t$とします.つまり,$X^{(1)},...,X^{(t)}$と$Y$(全区間)が観測されているとします.
  • 時刻$t$においてDTW距離が最も小さい区間の終了位置を$i^*$,DTW距離を$d_{min}$とします.この区間を候補と呼びます.

このとき,現在の観測時刻$t$において, 全ての$j = 1,...,T_y$について,$C_{t, j} \geq d_{min}$ または,$S_{t, j} \geq i^*$ が成り立つなら,現在の候補を検出結果として出力します.

これを図にすると以下です(突然の手書き図). 今,$t=5$で,$C_{1,T_y}$から$C_{t,T_y}$をチェックして,候補(赤線)が見つかったとします. このとき,SPRINGでは現在の時刻($t$)を通るパス(青や緑)について考えます. ここで,$t=6$以降はまだ観測していないので,これらのパスもDTW距離も計算することができないことに注意してください.

まず,緑色のパスは赤色の候補と区間が重複していないため,候補を検出するか否かを考える際に,考慮する必要がありません. $S_{t, j} \geq i^*$を満たすパスがこのパスに相当します.

次に青色のパスについて考えます. 青色のパスは候補と区間が重複しているため,もし青色のパスのほうがDTW距離が小さい場合は赤色の候補を検出してはいけません. これを判定する基準として,SPRINGでは$C_{t, j}$を用います. DTWでは正のフレーム間コストを考えるので,そのパスが高さ方向$j$の時点で$C_{t, j}$が$d_{min}$を超えた場合は,現在の候補よりも青色のパスのほうが必ずDTW距離が大きくなることがわかるからです.

f:id:ksknw:20191228103403j:plain

実際にSPRINGによって複数の区間を検出した結果を以下に示します.

def spring(x, y, epsilon):
    Tx = len(x)
    Ty = len(y)

    C = np.zeros((Tx, Ty))
    B = np.zeros((Tx, Ty, 2), int)
    S = np.zeros((Tx, Ty), int)

    C[0, 0] = dist(x[0], y[0])

    for j in range(1, Ty):
        C[0, j] = C[0, j - 1] + dist(x[0], y[j])
        B[0, j] = [0, j - 1]
        S[0, j] = S[0, j - 1]
        
    for i in range(1, Tx):
        C[i, 0] = dist(x[i], y[0])
        B[i, 0] = [0, 0]
        S[i, 0] = i
        
        for j in range(1, Ty):
            pi, pj, m = get_min(C[i - 1, j],
                                C[i, j - 1],
                                C[i - 1, j - 1],
                                i, j)
            C[i, j] = dist(x[i], y[j]) + m
            B[i, j] = [pi, pj]
            S[i, j] = S[pi, pj]
            
        imin = np.argmin(C[:(i+1), -1])
        dmin = C[imin, -1]
        
        if dmin > epsilon:
            continue
            
        for j in range(1, Ty):
            if (C[i,j] < dmin) and (S[i, j] < imin):
                break
        else:
            path = [[imin, Ty - 1]]
            temp_i = imin
            temp_j = Ty - 1
            
            while (B[temp_i, temp_j][0] != 0 or B[temp_i, temp_j][1] != 0):
                path.append(B[temp_i, temp_j])
                temp_i, temp_j = B[temp_i, temp_j].astype(int)
                
            C[S[:,-1] <= imin, -1] = 100000000
            yield np.array(path), dmin
pathes = []
for path, cost in spring(X, Y, 80):
    for line in path:
        plt.plot(line, [X[line[0]], Y[line[1]]], linewidth=0.8, c="gray")
    plt.plot(X)
    plt.plot(Y)
    plt.plot(path[:,0], X[path[:,0]], C="C2")
    plt.show()
    pathes.append(path)

f:id:ksknw:20191228103137p:plain

f:id:ksknw:20191228103151p:plain

path, cost = partial_dtw(X, Y)
D = (np.array(X).reshape(1, -1) - np.array(Y).reshape(-1, 1))**2
plot_path(pathes, X, Y, D)

f:id:ksknw:20191228103212p:plain

まとめ

ここでは,時系列データから特定の時系列パターンに一致した区間を検出する問題について考え,これを解くアルゴリズムの1つであるSPRINGをpythonで実装しました. 部分一致DTWは応用先の広いアルゴリズムで,最近ではDTWNet(ニューラルネットワークconvolutionカーネルの代わりにDTWを使う)に使われていました. 今回はpythonで実装したので遅いです.せっかくのオンライン設定用のアルゴリズムなので,気が向いたらjitとかjuliaとか使って高速化したいです.

参考

Blender + Pythonでポイントクラウドを可視化する

はじめに

ポイントクラウドデータをいい感じに可視化したい. matplotlibでも3次元データのscatterを描くことができるが,以下のような感じでいまいちな見た目になってしまう.

f:id:ksknw:20191029190640p:plain

もうちょっといい感じの図が描きたい.たとえはこんな感じのやつ. 調べてみると,blenderpythonが使えるらしいので,blenderなんもわからんけどやってみる. データはここで使われているShapenetの一部を用いた.

環境

ちょっと調べてでてくるスクリプトは動かないことが多かった(特にライト周り). 今使っているバージョン付近で大きめのアップデートがあったのか,blender自体がそういうアップデートの方針なのかわからないが,使っているバージョンによっては以下のスクリプトは動かないので注意.

Blenderのインストール

blenderはaptでもインストールできるが,最新版ではなかったのでsnapdというやつを使ってみることにした. 一般にパッケージ管理ソフトを複数混ぜるのは良くない気がするので,あまり良くない方法かもしれない.

sudo apt install snapd 
sudo snap install blender

terminalから

blender

で起動できる.

スクリプト

Blenderの中のScriptingタブからテキストエディタを開くことができる. Blenderから実行するpythonからはbpyというモジュールを使うことができ,これを使ってblenderの中のオブジェクトやカメラなどを操作できる. 全体はここに置いた.

オブジェクトの作成

ポイントクラウドを可視化するために,ポイントクラウドの各点に球を配置する. これは以下の手順でできる.

  1. 各位置にプリミティブを配置する.
  2. プリミティブの大きさを変更する.
  3. マテリアルを設定し,色を指定する.

プリミティブにはcubeやconeなど色々あるが,今回はuv_sphereを使う.

bpy.ops.mesh.primitive_uv_sphere_add(location=(x, y, z))

大きさの変更はオブジェクトを選択してスケールを変更することでできる. このあたりはかなりGUIpythonで操作している感じがある. bpy.context.object で作成したばかりのオブジェクト(選択済みになっているもの)を取得できる. 選択したオブジェクトのスケールを変更することで,球の大きさを変更できる.

sphere = bpy.context.object
sphere.scale = (0.02, 0.02, 0.02)

最後にマテリアルを追加することで,球の色を変更する. この機能を使うと,金属にするとか,水滴にするとかができる(と思う)が,今回は色を変更するだけ.

mat = bpy.data.materials.new("BLUE")
mat.diffuse_color = (.33, .43, .95, 1) # RGBA
sphere.data.materials.append(mat)

照明の配置

前述の通り、検索して出てくる多くのスクリプトが動かなかった。 これを参考にして,4方位から照らすようにした.

light_data = bpy.data.lights.new(name="light", type='AREA')
light_data.energy = 500
light_object = bpy.data.objects.new(name="light", object_data=light_data)
bpy.context.collection.objects.link(light_object)
bpy.context.view_layer.objects.active = light_object
light_object.location = loc
light_object.rotation_euler = rot

その他

カメラの配置はblender上で頑張って操作して決めた. renderのオプションなどは呪文として書いた.よくわかってない. 雰囲気的に床があったほうが良さそうだったので,適当に巨大な板を下に置いた.

bpy.ops.mesh.primitive_plane_add(location=(0, 0, -1))
plane = bpy.context.object
plane.scale = (100, 100, 1)

ということで以下のようなものを書いた.

実行

Blenderの中から実行しても良いし,コマンドラインから実行してもよい. jupyter notebookの中からも実行できるらしいが,まだやっていない (おそらく/snap/blender/33/2.80/python以下にipythonを入れればできそう).

コマンドラインから実行する場合は以下.

blender --python visualize_points_blender.py

--backgroundオプションをつけるとbackgroundで実行できるらしいが,セグフォでうまくできなかった.

実行すると以下のような図を得る.

f:id:ksknw:20191029190609p:plain

まとめ

Blenderpythonから使ってポイントクラウドを可視化した. GPUが弱いからか,1000点ぐらいのポイントクラウドを可視化するのに3分ぐらいかかっている. いい感じの図を作るセンスがほしい.

参考

Dockerを使ってpytorch+cudaの環境を構築する

はじめに

研究をやっていると

  • 査読対応のために実験を追加したいが,環境をいじってしまっている
  • 既存手法を実行するために,tensorflowの特定のバージョンを使いたい.それに伴って特定のバージョンのcudaをインストールしなければならない

みたいなことが発生する. pythonだけならcondaやpipenvでなんとかなるが,cudaなどが絡んでくると色々めんどくさくなる.

そこで,Dockerを使って環境を構築することでこの問題を解決することを試みる. もう1台のパソコンや会社のPCで同じことをやるので記録を残す.

cuda+dockerで検索するとnvidia-docker2を使うように書いてある記事がたくさんヒットするが,公式によると,nvidia-docker2はdepercatedでDocker 19.03以降はnvidia-container-toolkitを入れればいいらしい.

環境

家のパソコン

  • Ubuntu18.04
  • GeForce GTX 750 Ti, Core(TM) i7-7700 CPU
  • Docker version 19.03.1, build 74b1e89

gpuが年季入ってる感じなので不安だったけど大丈夫だった.

作りたい環境

  • cuda+pytorchの環境をjupyter notebook上で実行する
  • 作ったファイルはホスト上から参照,編集できるようにする

以下では,ほぼクリーンなubuntuに順番に環境を作っていく.

nvidia-driverのインストール

Dockerを使う場合でも,ホスト側にnvidiaのドライバを入れておく必要がある. 今どきはaptでいれると普通に動いてくれるので良い. ubuntu-driversでrecommendされたバージョンをインストールする.

sudo add-apt-repository ppa:graphics-drivers/ppa
sudo apt update
ubuntu-drivers devices

sudo apt install nvidia-driver-430 # recommended

参考:

Dockerのインストール

公式Documentに従ってインストールしていく.

sudo apt install apt-transport-https ca-certificates curl gnupg-agent software-properties-common
curl -fsSL https://download.docker.com/linux/ubuntu/gpg | sudo apt-key add -
sudo add-apt-repository "deb [arch=amd64] https://download.docker.com/linux/ubuntu $(lsb_release -cs) stable"
sudo apt install docker-ce docker-ce-cli containerd.io

sudo docker info

参考:

docker groupの設定

docker groupにuserを追加しておくと,管理者権限なしにdockerを実行できるようになる. dockerグループにはセキュリティ上の問題があると言っている記事もあるので,やらないほうがいいのかもしれないが理解できていないので,あまり検討できてない.

sudo gpasswd -a $USER docker
sudo reboot 

本当はログアウトして再度ログインするだけでいいはずだけど,自分の環境では再起動するまでgroupがうまく変更されていなかった.

docker hubのアカウント作成/ログイン

docker hubには様々なレポジトリが公開されており,これを利用することで誰かが作った色々盛り盛りの環境をさくっと利用することができる. 今回はnvidiaが公開しているものが使いたいので,これに登録,ログインしておく. ログインは以下のようにしてコマンドラインから実行する.

docker login

nvidia-container-toolkit のインストール

cuda+dockerで検索するとnvidia-docker2を使うように書いてある記事がたくさんヒットするが,公式によると,nvidia-docker2はdepercatedでDocker 19.03以降はnvidia-container-toolkitを入れればいいらしい.

Note that with the release of Docker 19.03, usage of nvidia-docker2 packages are deprecated since NVIDIA GPUs are now natively supported as devices in the Docker runtime.

Documentに従って,nvidia-container-toolkitをインストールする.

distribution=$(. /etc/os-release;echo $ID$VERSION_ID)
curl -s -L https://nvidia.github.io/nvidia-docker/gpgkey | sudo apt-key add -
curl -s -L https://nvidia.github.io/nvidia-docker/$distribution/nvidia-docker.list | sudo tee /etc/apt/sources.list.d/nvidia-docker.list

sudo apt-get update && sudo apt-get install -y nvidia-container-toolkit
sudo systemctl restart docker

docker build/run

以下では色々書いているけどnvidiaが公開しているDockerイメージを使うと万事オッケーかもしれない

Dockerfileを以下のように作る.

FROM nvidia/cuda:10.1-cudnn7-devel-ubuntu18.04

ARG USERNAME

RUN apt-get update
RUN apt-get -y upgrade
RUN apt-get -y install python3
RUN apt-get -y install python3-pip
RUN apt-get -y install nano wget curl

RUN pip3 install jupyter click numpy matplotlib seaborn pandas tqdm
RUN pip3 install torch torchvision

RUN useradd -m -s /bin/bash ${USERNAME}

USER ${USERNAME}

RUN jupyter notebook --generate-config \
 && sed -i.back \
    -e "s:^#c.NotebookApp.token = .*$:c.NotebookApp.token = u'':" \
    -e "s:^#c.NotebookApp.ip = .*$:c.NotebookApp.ip = '0.0.0.0':" \
    -e "s:^#c.NotebookApp.open_browser = .*$:c.NotebookApp.open_browser = False:" \
    /home/${USERNAME}/.jupyter/jupyter_notebook_config.py

WORKDIR /home/${USERNAME}

参考:

これをおいたフォルダで以下のようにしてイメージをビルドする.

docker build ./ --build-arg USERNAME=$USERNAME  -t notebook

ここではbuildの際にホスト側と同じユーザー名のユーザーを作っている. これによって,docker内で作ったファイルが,ホスト側で編集できるようになる. こちらにやると,別のマシンでビルドしたものを使うときにはこれだとだめらしいが,自分のマシンでしか使わないので,これでオッケー -tはイメージの名前なので何でも良い.

RUNしてみて,コンテナからGPUが見えることを確認する.

docker run --gpus all --rm notebook nvidia-smi

f:id:ksknw:20190831211946p:plain

ちゃんと動いていることが確認できた.

次にbashを実行してみる.

  • -itオプションをつけると,docker上でインタラクティブに操作できる.
  • また,-vでホームディレクトリをdocker内のhostディレクトリにマウントしておく.
  • -pはdockerのポートをホストのポートに飛ばすオプションで,jupyterを起動するポートをホストに飛ばしておく. これによって,docker上で起動したjupyterにホスト側からアクセスできるようになる.
docker run --gpus all --rm -v /home/kei:/home/kei/host -p 127.0.0.1:8888:8888 -it notebook bash

あとは普通にjupyter notebookを立ち上げると,いつもと同じようにホスト側のブラウザからlocalhost:8888にアクセスするだけで操作できる. bashを実行しているので,pythonを実行してもよい.

一応動作確認.

f:id:ksknw:20190831215419p:plain

家でやる場合はないけど,リモートサーバーでたてたdocker上のjupyter notebookにアクセスするときは,sshで入るときにポートフォワードすればよい.

参考:

おわりに

  • dockerを使うのはほぼ初めてなので,何か間違っているかもしれない.
  • 毎回--rmをつけてコンテナを使い捨てにするのが再現性的にはいいかなと思っている.
  • nvidia-docker2周りが最近変わったばかりっぽく,nvidia-docker2を入れろと書いてあるブログが多々あって混乱した.
  • これからはプロジェクトごとにDockerfileを作る感じでやっていきたい.

参考

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):な感じで入れておけば,必要な変数だけとってきて,残りはおいておくことができるので,いい感じに書けた.