1.はじめに

この記事はKerasのLSTMのフィードフォワードをnumpyで実装するの続きみたいなものです．
KerasでLSTM AutoEncoderを実装し，得られた特徴量から2値分類を試します．
データは，周波数の異なる2つのsin波を生成し，それを識別します．

2.環境

Python 3.6.1 :: Anaconda 4.4.0 (64-bit)
Keras 2.0.5 (backend : tensorflow)

3.データの作成

3.1 概要

データの生成は以下の記事を参考にしています．ありとうございます．

今回は2値分類をするため，データの生成方法を少し変更しています．
周波数の異なるsin波を二つ準備し，それを識別させます．

3.2 コード

今回使うデータを生成するコードです．

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

def make_data(random_factor, number_of_cycles, \
                timesteps, sampling_num_pair):
    """
    sampling_num_pair : 2 elements of tupple
        1周期の間にサンプリングする数
        Ex. (20,80)
    """
    def _load_data(data, n_prev = 100):  
        """
        data should be pd.DataFrame()
        """
        docX, docY = [], []
        for i in range(len(data)-n_prev):
            docX.append(data.iloc[i:i+n_prev].as_matrix())
        alsX = np.array(docX)

        return alsX

    def _train_test_split(df, test_size=0.1, n_prev = 100):  
        """
        This just splits data to training and testing parts
        """
        ntrn = round(len(df) * (1 - test_size))
        ntrn = int(ntrn)
        X_train = _load_data(df.iloc[0:ntrn], n_prev)
        X_test = _load_data(df.iloc[ntrn:], n_prev)

        return X_train, X_test

    np.random.seed(0)

    sampling_num1, sampling_num2 = sampling_num_pair

    df1 = pd.DataFrame(np.arange(sampling_num1 * number_of_cycles + 1), columns=["t"])
    df1["sin_t"] = df1.t.apply(lambda x: np.sin(x * (2 * np.pi / sampling_num1)+ np.random.uniform(-1.0, +1.0) * random_factor))

    df2 = pd.DataFrame(np.arange(sampling_num2 * number_of_cycles + 1), columns=["t"])
    df2["sin_t2"] = df2.t.apply(lambda x: np.sin(x * (2 * np.pi / sampling_num2)+ np.random.uniform(-1.0, +1.0) * random_factor))

    X_train1, X_test1 = _train_test_split(df1[["sin_t"]], n_prev=timesteps) 
    X_train2, X_test2 = _train_test_split(df2[["sin_t2"]], n_prev=timesteps) 

    # concatenate　X and make y
    X_train = np.r_[X_train1, X_train2]
    y_train = np.r_[np.tile(0, X_train1.shape[0]), np.tile(1, X_train2.shape[0])]

    X_test = np.r_[X_test1, X_test2]
    y_test = np.r_[np.tile(0, X_test1.shape[0]), np.tile(1, X_test2.shape[0])]

    return X_train, y_train, X_test, y_test

# 乱数の係数
random_factor = 0.05
# 生成するサイクル数
number_of_cycles = 200
# 1サイクルの中で，何点サンプリングするか．
sampling_num_pair=(20,80)
# windowの長さ．一つの系列の長さになる．
timesteps = 100

X_train, y_train, X_test, y_test = make_data(random_factor, number_of_cycles, timesteps, sampling_num_pair)
print("X_train.shape : ", X_train.shape) #X_train.shape :  (17802, 100, 1)
print("y_train.shape : ", y_train.shape) #y_train.shape :  (17802,)
print("X_test.shape : ", X_test.shape) #X_test.shape :  (1800, 100, 1)
print("y_test.shape : ", y_test.shape) #y_test.shape :  (1800,)

3.3 プロット

生成したsin波をプロットしてみます．

# make_data関数のdf1, df2を抜き出してプロットする．
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(df1["sin_t"][0:80], label="class0 sampling_num 20", color="red")
ax.plot(df2["sin_t2"][0:80], label="class1 sampling_num 80", color="blue")
ax.legend(loc="upper right")
fig.savefig("./sin_plot.png")
plt.close("all")

class0はclass1の4倍の周波数であることがわかります．

4.model構築

4.1 概要

モデルの構築にはThe Keras Blogを参考にしました．Sequence-to-sequence autoencoderの章にヒントっぽいものがあったのでそれに適宜追加しています．

元になる論文はUnsupervised Learning of Video Representations using LSTMsです．この論文ではcondtional（decode時に入力を逆順にしたものを渡す）とunconditinal（何も渡さない）を提案していますが，今回はcondtionalの方を実装しています．

4.2 コード

LSTM AutoEncoderのモデルを定義します．
ついでに訓練させて，元のX_trainを再構成します．

from keras.layers import Input, LSTM, RepeatVector, concatenate, Dense
from keras.models import Model

input_dim = 1
latent_dim = 10

# encode
inputs = Input(shape=(timesteps, input_dim))
encoded = LSTM(latent_dim, activation="tanh", recurrent_activation="sigmoid", return_sequences=False)(inputs)

#decode
hidden = RepeatVector(timesteps)(encoded)
reverse_input = Input(shape=(timesteps, input_dim))
hidden_revinput = concatenate([hidden, reverse_input])
decoded = LSTM(latent_dim, activation="tanh", recurrent_activation="sigmoid", return_sequences=True)(hidden_revinput)
decoded = Dense(latent_dim, activation="relu")(decoded)
decoded = Dense(input_dim, activation="tanh")(decoded)

# train
LSTM_AE = Model([inputs, reverse_input], decoded)
LSTM_AE.compile(optimizer='rmsprop', loss='mse')
X_train_rev = X_train[:,::-1,:]
LSTM_AE.fit([X_train, X_train_rev], X_train, epochs=30, batch_size=500, shuffle=True, validation_data=([X_train, X_train_rev], X_train))

X_hat = LSTM_AE.predict([X_train, X_train_rev])

このコードではdecode時にフルコネクト層を追加していますが，特に深い意味はありません．Kerasすげえ！こんな簡単に層を追加できるんだ！と遊んでいた結果です．

次に，class0とclass1について，どのようにAutoEncoderで再構成されているか確認するため，classごとに分割します．

def split_X(X, y):
    y_inv = np.abs(y - 1.)
    X_0 = X[y_inv.astype(bool),:,:]
    X_1 = X[y.astype(bool),:,:]
    return X_0, X_1

X_train_0, X_train_1 = split_X(X_train, y_train)
X_hat_0, X_hat_1 = split_X(X_hat, y_train)

print("X_train_0.shape : ", X_train_0.shape) #X_train_0.shape :  (3501, 100, 1)
print("X_train_1.shape : ", X_train_1.shape) #X_train_1.shape :  (14301, 100, 1)
print("X_hat_0.shape : ", X_hat_0.shape) #X_hat_0.shape :  (3501, 100, 1)
print("X_hat_1.shape : ", X_hat_1.shape) #X_hat_1.shape :  (14301, 100, 1)

4.3 プロット

AutoEncoderで再構成されたものをプロットして確認します．

# reconstruct したX_trainがどんな感じか見てみる
def plot_save(start_index, X_hat, X_train, X_class):
    fig = plt.figure()
    ax = fig.add_subplot(111)
    for i in np.arange(start_index,start_index+5):
        #５個ずつプロットする．
        ax.plot(X_hat[i,:,0], label="X hat", color="red")
        ax.plot(X_train[i,:,0], label="X train", color="blue")
    savename = "./AE_reconst_30ep_start_" + str(start_index) + "_cls" + str(X_class) + ".png"
    fig.savefig(savename)
    plt.close("all")

start_list = np.arange(0,len(X_train_0), 1000)
for start_index in start_list:
    plot_save(start_index, X_hat_0, X_train_0, X_class=0)

start_list = np.arange(0,len(X_train_1), 1000)
for start_index in start_list:
    plot_save(start_index, X_hat_1, X_train_1, X_class=1)

class0

class1

青が元のデータ，赤が再構築したものです．
これらは，何枚もプロットした中で適当に選んだものです．どれもだいたいこんな感じでした．

ちなみにまだ収束していない様子だったので，epochを増やせばもっと良い近似になると思いますが，マシンスペック的にまあこのくらいでいいだろうということにしました．ちなみに30epochで15分ほどでした．（もっと頑張れよと言われそう）

5.獲得した特徴量の抽出

5.1 概要

このへんは前回の記事と同じです．
encoderで得た特徴量をmodelから取り出します．

5.2 コード

特徴量を抜き出すのに必要な関数を定義します．


def split_params(W, U, b):
    Wi = W[:,0:latent_dim]
    Wf = W[:,latent_dim:2*latent_dim]
    Wc = W[:,2*latent_dim:3*latent_dim]
    Wo = W[:,3*latent_dim:]

    print("Wi : ",Wi.shape)
    print("Wf : ",Wf.shape)
    print("Wc : ",Wc.shape)
    print("Wo : ",Wo.shape)

    Ui = U[:,0:latent_dim]
    Uf = U[:,latent_dim:2*latent_dim]
    Uc = U[:,2*latent_dim:3*latent_dim]
    Uo = U[:,3*latent_dim:]

    print("Ui : ",Ui.shape)
    print("Uf : ",Uf.shape)
    print("Uc : ",Uc.shape)
    print("Uo : ",Uo.shape)

    bi = b[0:latent_dim]
    bf = b[latent_dim:2*latent_dim]
    bc = b[2*latent_dim:3*latent_dim]
    bo = b[3*latent_dim:]
    print("bi : ",bi.shape)
    print("bf : ",bf.shape)
    print("bc : ",bc.shape)
    print("bo : ",bo.shape)

    return (Wi, Wf, Wc, Wo), (Ui, Uf, Uc, Uo), (bi, bf, bc, bo)


def calc_ht(params):
    x, latent_dim, W_, U_, b_ = params
    Wi, Wf, Wc, Wo = W_
    Ui, Uf, Uc, Uo = U_
    bi, bf, bc, bo = b_ 
    n = x.shape[0]

    ht_1 = np.zeros(n*latent_dim).reshape(n,latent_dim) #h_{t-1}を意味する．
    Ct_1 = np.zeros(n*latent_dim).reshape(n,latent_dim) #C_{t-1}を意味する．

    ht_list = []

    for t in np.arange(x.shape[1]):
        xt = np.array(x[:,t,:])
        it = sigmoid(np.dot(xt, Wi) + np.dot(ht_1, Ui) + bi)
        Ct_tilda = np.tanh(np.dot(xt, Wc) + np.dot(ht_1, Uc) + bc)
        ft = sigmoid(np.dot(xt, Wf) + np.dot(ht_1, Uf) + bf)
        Ct = it * Ct_tilda + ft * Ct_1
        ot = sigmoid( np.dot(xt, Wo) + np.dot(ht_1, Uo) + bo)
        ht = ot * np.tanh(Ct)
        ht_list.append(ht)
        ht_1 = ht
        Ct_1 = Ct

    ht = np.array(ht)
    return ht
    
def sigmoid(x):
    return 1.0 / (1.0 + np.exp(-x))

実際に抜き出します．

W, U, b = LSTM_AE.layers[1].get_weights()
Ws, Us, bs = split_params(W, U, b)

params = [X_train, latent_dim, Ws, Us, bs]
ht_train = calc_ht(params)

params = [X_test, latent_dim, Ws, Us, bs]
ht_test = calc_ht(params)

LSTM_AE.layers[1]はencodeのLSTM層です．そこから重みを抜き出しています．

5.3プロット

ht_trainは (17802, 10) の行列で，10次元の特徴量になっています.
元の次元がtimesteps = 100で設定した100次元なので，100次元から10次元に圧縮されています．
この10次元の特徴量をとりあえずプロットしてみましょう．

class0

class1

なんだかよく分からないですね．周期的な信号なので，1周期の特徴量の動きが全て表示されるようにプロットしています．全体的な挙動としてはだいたい似ているけれども，まあ一個ずつみたら違うかもなあと言う感じですね．（大変なので細かくはみない）

6.獲得した特徴量で識別

このLSTM AutoEncoderによって獲得した特徴量によってclass0とclass1を分類してみましょう．とりあえずここでは，iidじゃないじゃんと言うことは置いておいて，ロジスティック回帰を試します．

from sklearn.linear_model import LogisticRegression 
from sklearn.metrics import accuracy_score
model = LogisticRegression(n_jobs=-1)
model.fit(ht_train, y_train)
y_hat_test = model.predict(ht_test)
accuracy_score(y_hat_test, y_test)
# 1.0

テストデータに対するaccuracyが1.0になりました．高すぎて逆に心配なのですが，かなり波形の違う信号でしたしこんなものなんでしょうか．

7.終わりに

Kerasで初めてモデル構築をしたが，高レイヤーすぎて中で実際にやっていることが自分のやりたいこととあっているのか心配になる．結果だけ見ると期待通りになっているので，できていると信じたい．（誰か教えてください）
これの元の論文もそうだが，そこで提案されているモデルが数式で書かれていなかったりして，機械学習界隈の人はどうやって理解したり実装したりしているのだろうか．computational graphという考え方が広がっていて，図によって伝達されているのだろうか．それとも，論文読んだだけではモデルが想像できず，再現性がないと問題になっているのだろうか．よく分からない．
Kerasはtendorflowのラッパーなので当然なのだが，tensolflow触ったことがない僕のような人間は困難にぶち当たるのでとりあえずtensorflowを触るのがいいかもしれない．