Kunst Site
Blogブログ
About Site当サイトについて
Portfolioポートフォリオ
Contact問い合わせ

GPLVMを実装してみた


データサイエンス教師なし学習次元削減
更新日:2022-02-22

GPLVM を python で実装をしてみました。

概要

ガウス過程潜在変数モデル(GPLVM)は、ガウス過程を使って(潜在的に)高次元データの低次元表現を教師なし学習する次元削減方法です。 **本記事では、ガウス過程の青本を参考にしています。

問題設定

NNN個の観測データをX={xn}n=1N,xn∈X\mathbf X = \{\mathbf x_n \}^N_{n=1}, \mathbf x_n \in \mathcal XX={xn​}n=1N​,xn​∈Xが与えられたとき、X\mathcal{X}X は高次元の観測空間で、DDD次元のユークリッド空間とします。つまりX=RD\mathcal X = \mathbb R^DX=RDであり、データ全体はX∈RN×D\mathbf X \in \mathbb R^{N \times D}X∈RN×Dと表せます。

X\mathbf{X}Xは各次元で別々の意味を持ち、近いz\mathbf{z}z同士では近い値になっていると考えます。ここで、NNN個のxn\mathbf{x_n}xn​のそれぞれに対応した未知のNNN個の潜在変数Z={zn}n=1N\mathbf{Z}=\{\mathbf{z_n}\}^N_{n=1}Z={zn​}n=1N​からのガウス過程回帰での写像fffによって生成されると仮定します。GPLVM は、この写像fffを推定します。

gplvm問題設定

実装

GPLVM では観測データ同士が近いところにあれば、それらに対応する潜在変数も近いところにあると考えます。 NNN個のxn\mathbf{x_n}xn​のそれぞれに対応した未知のNNN個の潜在変数{zn}n=1N\{\mathbf{z_n}\}^N_{n=1}{zn​}n=1N​からのガウス過程回帰での写像fffによって生成されると仮定しているので、Z\mathbf{Z}Zがわかれば、それぞれの各出力次元は独立なので、データXXX全体の確率は、各データの次元の和になります。

p(X∣Z)=∏d=1Dp(x∣Z)=∏d=1DN(x∣0,Kz+σ2I)p(\mathbf{X}|\mathbf{Z})=\prod_{d=1}^D p(\mathbf{x}|\mathbf{Z})=\prod_{d=1}^D N(\mathbf{x}|0, \mathbf{K_z}+\sigma^2 \mathbf{I})p(X∣Z)=d=1∏D​p(x∣Z)=d=1∏D​N(x∣0,Kz​+σ2I)

ここでx\mathbf{x}xは平均が000で、ガウス分布による誤差を考えて以下の式で表せられるとしています。また、それぞれのz\mathbf{z}zがKKK次元の標準ガウス分布に従うと考えます。

x∼(0,Kz+σ2I)\mathbf{x} \sim (0,\mathbf{K_z}+\sigma^2\mathbf{I})x∼(0,Kz​+σ2I) p(X)=∏n=1Np(zn),    p(x)=N(0,IK)p(\mathbf{X})=\prod_{n=1}^N p(\mathbf{z_n}),\,\,\,\,p(\mathbf{x})=N(0,\mathbf{I_K})p(X)=n=1∏N​p(zn​),p(x)=N(0,IK​)

これらのことからX\mathbf{X}XとZ\mathbf{Z}Zの同時確率は次のように表すことができます。

p(X,Z)=p(X∣Z)p(Z)=∏d=1DN(x∣0,Kz+σ2I)∏n=1NN(zn∣0,IK)p(\mathbf{X},\mathbf{Z})=p(\mathbf{X}|\mathbf{Z})p(\mathbf{Z})=\prod_{d=1}^DN(\mathbf{x}|0,\mathbf{K_z}+\sigma^2\mathbf{I})\prod_{n=1}^NN(\mathbf{z_n}|0,\mathbf{I_K})p(X,Z)=p(X∣Z)p(Z)=d=1∏D​N(x∣0,Kz​+σ2I)n=1∏N​N(zn​∣0,IK​)

GPLVM では、この式が最大となるようなZ\mathbf{Z}Zを見つけます。

詳細な計算は省きますが、この式を変形して対数尤度を考えると下記の式になります。このとき、青本では∏n=1NN(zn∣0,IK)\prod_{n=1}^NN(\mathbf{z_n}|0,\mathbf{I_K})∏n=1N​N(zn​∣0,IK​)の対数尤度(おそらく正則化項)は無視しています・

L=log⁡p(X∣Z)=−ND2log⁡(2π)−D2log⁡∣Kz∣−12tr(Kz−1YYT)L=\log p(\mathbf{X}|\mathbf{Z})=-\frac{ND}{2}\log(2\pi)-\frac{D}{2}\log|\mathbf{K_z}|-\frac{1}{2}tr(\mathbf{K_z^{-1}}\mathbf{Y}\mathbf{Y}^T)L=logp(X∣Z)=−2ND​log(2π)−2D​log∣Kz​∣−21​tr(Kz−1​YYT)

観測変数の対数尤度LLLから潜在変数の更新の勾配法を求められます。計算していくと下記の式が導けます。

∂L∂Kz=12(Kx−1YYTKz−1−DKz−1)\frac{\partial L}{\partial \mathbf{K_z}}=\frac{1}{2}(\mathbf{K^{-1}_x}\mathbf{Y}\mathbf{Y}^T\mathbf{K_z^{-1}}-D\mathbf{K_z^{-1}})∂Kz​∂L​=21​(Kx−1​YYTKz−1​−DKz−1​) ∂Kz∂x=−2k(zn,zn′)(znj−zn′j)/σ\frac{\partial \mathbf{K_z}}{\partial \mathbf{x}}=-2k(\mathbf{z_n},\mathbf{z}_n^{\prime})(z_{nj}-z_{n^{\prime}j})/\sigma∂x∂Kz​​=−2k(zn​,zn′​)(znj​−zn′j​)/σ

このとき、カーネル関数は下記のガウスカーネルを用いています。δ\deltaδのところはn=n′n=n^{\prime}n=n′ならば111となり、それ以外は000になります。

K(z,z′)=θ1exp⁡(−∣zn−zn′∣2θ2)+θ3δ(n,n′)K(z,z^{\prime})=\theta_1\exp\left(-\frac{|z_n-z_n^{\prime}|^2}{\theta_2}\right)+\theta_3\delta(n,n^{\prime})K(z,z′)=θ1​exp(−θ2​∣zn​−zn′​∣2​)+θ3​δ(n,n′)

潜在変数の更新と写像の推定

したがって∂L∂z=∂L∂Kz∂Kz∂x\frac{\partial L}{\partial \mathbf z}=\frac{\partial L}{\partial \mathbf{K_z}}\frac{\partial \mathbf{K_z}}{\partial \mathbf{x}}∂z∂L​=∂Kz​∂L​∂x∂Kz​​より、潜在変数の更新と写像の推定は以下の式で行います。

潜在変数の更新(勾配法)

∂L∂z\frac{\partial L}{\partial \mathbf z}∂z∂L​はマイナスなので下記の式になります。

z=z+∂L∂z\mathbf{z}=\mathbf{z}+\frac{\partial L}{\partial \mathbf z}z=z+∂z∂L​

写像の推定

GPR を用いています。青本では太字を横ベクトルで表していますが、実装するときは縦ベクトルで考えているので転置のところは消えています。

f(z)∼GP(K(z,z′)(K(z,z))−1X,    K(z′,z′)−K(z,z′)(K(z,z))−1K(z,z′))f(\mathbf{z}) \sim \mathcal{GP}(\mathbf{K}(\mathbf{z},\mathbf{z}^{\prime})(\mathbf{K}(\mathbf{z},\mathbf{z}))^{-1}\mathbf{X},\,\,\,\,\mathbf{K}(\mathbf{z}^{\prime},\mathbf{z}^{\prime})-\mathbf{K}(\mathbf{z},\mathbf{z}^{\prime})(\mathbf{K}(\mathbf{z},\mathbf{z}))^{-1}\mathbf{K}(\mathbf{z},\mathbf{z}^{\prime}))f(z)∼GP(K(z,z′)(K(z,z))−1X,K(z′,z′)−K(z,z′)(K(z,z))−1K(z,z′))

実装コードと結果

import numpy as np
import matplotlib.animation as animation
import matplotlib.pyplot as plt
from scipy.spatial import distance as dist
from tqdm import tqdm

class GPLVM(object):
def **init**(self, θ1, θ2, θ3):
self.θ1 = θ1
self.θ2 = θ2
self.θ3 = θ3

    def fit(self, X, latent_dim, epoch, eta):
        resolution = 10
        T = epoch
        N, D = X.shape
        L = latent_dim
        Z = np.random.randn(N, L) /100

        history = {}
        history['Z'] = np.zeros((T, N, L))
        history['f'] = np.zeros((T, resolution, resolution, D))

        for t in tqdm(range(T)):
            K = self.θ1 * self.kernel(Z, Z, self.θ2) +  self.θ3 * np.eye(N)
            inv_K = np.linalg.inv(K)
            dLdK = 0.5 * (inv_K @ (X @ X.T) @ inv_K - D * inv_K)
            dKdX = -2.0*(((Z[:,None,:]-Z[None,:,:])*K[:,:,None]))/self.θ2
            dLdX = np.sum(dLdK[:, :, None] *  dKdX, axis=1)

            Z = Z + eta * dLdX
            history['Z'][t] = Z

            z_new_x = np.linspace(min(Z[:,0]),max(Z[:,0]), resolution)
            z_new_y = np.linspace(min(Z[:,1]),max(Z[:,1]), resolution)
            z_new = np.dstack(np.meshgrid(z_new_x, z_new_y)).reshape(resolution**2, L)
            k_star = self.θ1 * self.kernel(z_new, Z, self.θ2)
            F = (k_star @ inv_K @ X).reshape(resolution, resolution, D)
            history['f'][t] = F
        return history

    def kernel(self,X1, X2, θ2):
        Dist = np.sum(((X1[: , None, :] - X2[None, :, :])**2), axis=2)
        K = np.exp((-0.5/θ2) * Dist)
        return K

if **name** == "**main**":
np.random.seed(0)
resolution = 100
z1 = np.random.rand(resolution) _ 2.0 - 1.0
z2 = np.random.rand(resolution) _ 2.0 - 1.0
X = np.empty((resolution, 3))
X[:, 0] = z1
X[:, 1] = z2
X[:, 2] = (z1**2 - z2**2)
 X += np.random.normal(loc=0, scale=0.0, size=X.shape)

    model = GPLVM(θ1=1.0, θ2=0.03, θ3=0.05)
    history = model.fit(X,latent_dim=2, epoch=100, eta=0.00001)

# ---------描写---------------------------------------------------------------

    fig = plt.figure(figsize=(10, 5))
    ax_observable = fig.add_subplot(122, projection='3d')
    ax_latent = fig.add_subplot(121)

    def update(i, z, x, f):
        plt.cla()
        ax_latent.cla()
        ax_observable.cla()

        fig.suptitle(f"epoch: {i}")
        ax_latent.scatter(z[i,:, 0], z[i,:, 1],s=50, edgecolors="k", c=x[:,0])
        ax_observable.scatter(x[:, 0], x[:, 1], x[:, 2], c=x[:,0],s=50, marker='x')
        ax_observable.plot_wireframe(f[i ,:, :, 0], f[i, :, :, 1], f[i, :, :, 2], color='black')

        ax_observable.set_xlim(x[:, 0].min(), x[:, 0].max())
        ax_observable.set_ylim(x[:, 1].min(), x[:, 1].max())
        ax_observable.set_zlim(x[:, 2].min(), x[:, 2].max())

        ax_observable.set_title('observable_space')
        ax_observable.set_xlabel("X_dim")
        ax_observable.set_ylabel("Y_dim")
        ax_observable.set_zlabel("Z_dim")
        ax_latent.set_title('latent_space')
        ax_latent.set_xlabel("X_dim")
        ax_latent.set_ylabel("Y_dim")

    ani = animation.FuncAnimation(fig, update, fargs=(history['Z'], X, history['f']), interval=100, frames=100)
    # ani.save("tmp.gif", writer = "pillow")
    print("X: {}, Y: {}, Z:{}".format(X.shape, history['f'][0].shape, history['Z'][0].shape))
    plt.show()
result

右図が観測空間上での推定した多様体の学習過程で、右図はその時の潜在変数です。ハイパーパラメータを手動で探しましたが、綺麗に描画してくれるのを見つけるのに苦労しました。コンマ単位でハイパーパラメータを変えると描画ができなくなってしまいます。GPLVM がそれだけ繊細なモノということを実感しました。

参考文献

ガウス過程と機械学習 (機械学習プロフェッショナルシリーズ) 持橋 大地 (著), 大羽成征

終わりに

ここまで読んでくださってありがとうございます。

編集リクエストもお待ちしてます。

Profile

kunst

データサイエンティスト兼フロントエンジニアのkunstです!
最近は山登りにハマっています。

目次

kunst Site

Welcome to kunst Site. Stay updated with the latest posts.

Privacy PolicyAbout Blog|BlogAbout SitePortfolioContact

Copyright © 2024 kunst Site