pythonコードを修正しました．
この記事は古川研究室Advent_calendar 23日目の記事です。

はじめに

GPLVMはガウス過程(GP)を用いた教師なし学習法で，高次元空間上にデータが分布する低次元の多様体を推定することができます．GPLVMの一番の魅力はモデルがとてもシンプルであるため拡張性が高く，理論的に扱いやすいところです．実際，ハイパーパラメータをベイズ推定の枠組みで推定したり，時系列データ解析¹やマルチビュー解析²，テンソル分解³などのより複雑な状況への拡張されてきました．今回はもっとも基本的なGPLVMのアルゴリズムを確率的主成分分析(probabilistic PCA : pPCA)からスタートして導出します．今回の導出はおおむねLawrence先生の論文の流れに従っているので詳細はそちらの方をご参照ください⁴．

probabilistic PCA

pPCAはその名の通り主成分分析を確率論の枠組みで主成分分析を再導出したモデルです． pPCAでは$\mathbf{X}=(\mathbf{x}_1,\mathbf{x}_2,\cdots,\mathbf{x}_N)\in\mathbb{R}^{D \times N}$が与えられた時に潜在変数$\mathbf{z} \in \mathbb{R}^L$が以下のように写像されることで$\mathbf{x}$が観測されると仮定します⁵．

\mathbf{x} = \mathbf{W}\mathbf{z} + \boldsymbol{\epsilon},

ここで$\epsilon$は観測ノイズであり精度パラメータが$\beta^{-1}$の等方ガウスノイズを仮定します．また$\mathbf{W}\in \mathbb{R}^{D \times L}$です．また$\mathbf{z}$は$L$次元の標準正規分布に従うと仮定します．つまり$\mathbf{x}$と$\mathbf{z}$は次のような確率分布に従います．

\begin{align}
p(\mathbf{x}\mid\mathbf{W},\mathbf{z},\beta) &= \mathcal{N}(\mathbf{x}\mid\mathbf{W}\mathbf{z},\beta^{-1}\mathbf{I}_D) \\
p(\mathbf{z}) &= \mathcal{N}(\mathbf{z}\mid\mathbf{0},\mathbf{I}_L)
\end{align}

この時$p(\mathbf{x}\mid\mathbf{W},\mathbf{z},\beta)$も$p(\mathbf{z})$もどちらもガウス分布であるため同時分布$p(\mathbf{x},\mathbf{z}\mid\mathbf{W},\beta)$もガウス分布となり，$p(\mathbf{x},\mathbf{z}\mid\mathbf{W},\beta)$を潜在変数$\mathbf{z}$について周辺化した$p(\mathbf{x}\mid\mathbf{W},\beta)$もまたガウス分布となります．実際これは$p(\mathbf{x}\mid\mathbf{W},\beta)=\mathcal{N}(\mathbf{x}\mid\mathbb{E}[\mathbf{x}],\mathbb{V}[\mathbf{x}])$とすると以下のように計算できます．

\begin{align}
\mathbb{E}[\mathbf{x}]&=\mathbb{E}[\mathbf{W}\mathbf{z} + \boldsymbol{\epsilon}] \\
&=\mathbf{W}\mathbb{E}[\mathbf{z}] + \mathbb{E}[\boldsymbol{\epsilon}] \\
&=\mathbf{0} \\
\mathbb{V}[\mathbf{x}]&=\mathbb{E}[(\mathbf{W}\mathbf{z} + \boldsymbol{\epsilon})(\mathbf{W}\mathbf{z} + \boldsymbol{\epsilon})^{\rm T}] \\
&=\mathbf{W}\mathbb{E}[\mathbf{z}\mathbf{z}^{\rm T}]\mathbf{W}^{\rm T} + \mathbb{E}[\boldsymbol{\epsilon}\mathbf{z}^{\rm T}]\mathbf{W}^{\rm T} + \mathbb{E}[\boldsymbol{\epsilon}\boldsymbol{\epsilon}^{\rm T}]\\
&=\mathbf{W}\mathbf{W}^{\rm T} + \beta^{-1}\mathbf{I}_D \\
&=\mathbf{C}
\end{align}

上記の式変形では$\mathbf{W}$は確率変数でないため期待値の外に出すことができること，$\mathbf{z}$と$\boldsymbol{\epsilon}$が独立と仮定していることを用いています．

ここで$\boldsymbol{\epsilon}$も$\mathbf{z}$も独立同分布であるため$p(\mathbf{x}\mid\mathbf{W},\beta)$も$\mathbf{x}$について独立同分布となります．したがって

\begin{align}
p(\mathbf{X}\mid\mathbf{W},\beta)=\prod^N_{n=1}p(\mathbf{x}_n\mid\mathbf{W},\beta)
\tag{2}
\end{align}

となります．以下のように式(2)の尤度関数の対数を$\mathbf{W}$と$\beta$について最大化を行います．

\begin{align}
L_{\rm pPCA}&=\log{p(\mathbf{X}\mid\mathbf{W},\beta)} \\
&=\sum^N_{n=1}\log{p(\mathbf{x}_n\mid\mathbf{W},\beta)} \\
&=\sum^N_{n=1}\left(-\frac{D}{2}\log{2\pi}-\frac{1}{2}\log{|\mathbf{C}|}-\frac{1}{2}{\rm Tr}[\mathbf{C}^{-1}\mathbf{x}_n\mathbf{x}^{\rm T}_n]\right) \\
&=-\frac{ND}{2}\log{2\pi}-\frac{N}{2}\log{|\mathbf{C}|}-\frac{N}{2}{\rm Tr}[\mathbf{C}^{-1}\mathbf{V}] \tag{3}
\end{align}

これが確率的主成分分析の目的関数となります⁶．この目的関数を最大にする$\mathbf{W}$と$\beta$を推定する方法は$L$を微分して閉形式の形で求める方法とEMアルゴリズムを用いて求める方法がありますが，本題から外れるためここでは省略します．詳細はこちらの論文を参照ください⁷

Dual probabilistic PCA

pPCAでは$\mathbf z_n$を確率変数，$\mathbf{W}$をパラメータとしてみなして$\mathbf x_n$の確率モデルを考えました．
これに対しDual Probabilistic PCA (DPPCA)では$\mathbf{Z}=(\mathbf z_1,\mathbf z_2,\cdots,\mathbf z_N)\in\mathbb{R}^{L\times N}$をパラメータ，$\mathbf w_d$を確率変数とみなして
$\mathbf x_{:d}$の確率モデルを考えます．ここで$\mathbf x_{:d}=(x_{1d},x_{2d},\cdots,x_{Nd})^{\rm T}$です．つまり，

\mathbf{x}_{:d} = \mathbf{Z}^{\rm T}\mathbf{w}_d + \boldsymbol{\epsilon}_{:d},

とし，$\mathbf w_d$と$\boldsymbol \epsilon_{:d}$はそれぞれ以下の確率分布に従うと仮定します．

\begin{align}
p(\mathbf w_d)&=\mathcal{N}(\mathbf w_d\mid\mathbf{0},\mathbf{I}_L) \\
p(\boldsymbol{\epsilon}_{:d})&=\mathcal{N}(\boldsymbol{\epsilon}_{:d}\mid\mathbf{0},\beta^{-1}\mathbf{I}_N)
\end{align}

$p(\mathbf w_d)$も$p(\mathbf x_{:d}\mid \mathbf{Z},\mathbf{w}_{:d},\beta)$もガウス分布なのであとはpPCAと同様に$\mathbf{w}_d$に関して周辺化することで尤度関数が以下のように求まります．

\begin{align}
p(\mathbf{X} \mid \mathbf{Z}, \beta)&=\prod^D_{d=1}\mathcal{N}(\mathbf{x}_{:d}\mid\mathbf{0},\mathbf{Z}\mathbf{Z}^{\rm T}+\beta^{-1}\mathbf{I}_N) \\
\end{align}

この尤度関数の対数をとるとDPPCAの目的関数が求まります．

\begin{align}
\log{p(\mathbf{X} \mid \mathbf{Z}, \beta)} &=\sum^D_{d=1}\log{p(\mathbf{x}_{:d} \mid \mathbf{Z}, \beta)} \\
&=-\frac{ND}{2}\log{2\pi}-\frac{D}{2}\log{|\mathbf{K}|}-\frac{D}{2}{\rm Tr}[\mathbf{K}^{-1}\mathbf{S}] \tag{4}
\end{align}

ここで$\mathbf S$と$\mathbf K$はそれぞれ以下のように定義しました．

\begin{align}
\mathbf S &= \frac 1 D \mathbf{X}\mathbf{X}^{\rm T} \\
\mathbf K &= \mathbf Z \mathbf Z^{\rm T}+\beta^{-1}\mathbf I_N 
\end{align}

式(4)を見ると結局のところ$\mathbf W$と$\mathbf Z$の役割が入れ替わっているだけなのでpPCAとやっていることはそれほど違いがないように思われるかもしれません．しかし以下の定数から式(4)を引いたものを考えてみます．

\int \mathcal{N}(\mathbf{x}\mid\mathbf{0},\mathbf{S})\log{\mathcal{N}(\mathbf{x}\mid\mathbf{0},\mathbf{S})}d\mathbf{x}=-\frac{ND}{2}\log{2\pi}-\frac{D}{2}\log{|\mathbf S|}+\frac{ND}{2}

この時，以下の式より式(4)を最大化する$\mathbf{Z}$を推定することは観測データのグラム行列と潜在変数のグラム行列の間のKLダイバージェンスを最小にする$\mathbf{Z}$を推定することと等価であることがわかります．つまり観測データ間の類似度と対応する潜在変数間の類似度がなるべく一致するように潜在変数$\mathbf{Z}$を推定することがDPPCAの目的と言えます．

\begin{align}
\int \mathcal{N}(\mathbf{x}\mid\mathbf{0},\mathbf{S})\log{\mathcal{N}(\mathbf{x}\mid\mathbf{0},\mathbf{S})}d\mathbf{x}-L_{\rm DPPCA} &=\frac{D}{2}\{-\log{|\mathbf S|}+\log{|\mathbf{K}|}+{\rm Tr}[\mathbf{K}^{-1}\mathbf{S}]+N\} \\
&= D_{\rm KL}[\mathcal{N}(\mathbf{x}\mid\mathbf{0},\mathbf{S}) || \mathcal{N}(\mathbf{x}\mid\mathbf{0},\mathbf{K})]
\end{align}

この時DPPCAでは観測データの類似度と潜在変数の類似度をカーネル関数$k(\cdot,\cdot)$を用いて標準内積以外の類似度を定義することで非線形な次元削減を実現できます．特に観測データの類似度をカーネル関数を用いて非線形な次元削減を行う方法をカーネル主成分分析と呼び⁸，潜在変数の類似度をカーネル関数を用いて非線形な次元削減を行う方法をGPLVMと呼びます．

それぞれの手法はどちらも非線形な次元削減が可能ですがそれぞれ長所と短所があります．カーネル主成分分析では観測データ間の類似度にカーネル関数を用いるためカーネル関数を用いて一度$\mathbf{S}$を計算してしまえば通常の主成分分析と同様に解析解が得られます．しかし，潜在空間から観測空間への写像がわからないため潜在空間上の任意の点に対応する観測空間上の点を推定するためにはpre-image問題⁹を解く必要が出てきます．これに対しGPLVMでは潜在空間から観測空間への写像を陽に書くことができるためpre-image問題は発生しません．その代わり潜在変数が変わるたびに$\mathbf{K}$を更新する必要があります．

Gaussian Process Latent Variable Model

これまでの話でGPLVMというのは以下の目的関数を最大にするように潜在変数$\mathbf Z$を推定する学習モデルであることがわかりました．

L_{\rm DPPCA} = -\frac{ND}{2}\log{2\pi}-\frac{D}{2}\log{|\mathbf{K}|}-\frac{D}{2}{\rm Tr}[\mathbf{K}^{-1}\mathbf{S}] \tag{5}

ここで$\mathbf K$はカーネル関数$k(\cdot,\cdot)$を用いて$\mathbf K= k(\mathbf Z,\mathbf Z) + \beta^{-1}\mathbf{I}_N$と定義されます．

実はこの目的関数はGPの観点で導出することもできます．今回はGPについては詳細は省略しますのでGPについて知りたい方はこちらの本をご参照ください¹⁰．潜在空間$\mathcal Z$から観測空間$\mathcal X$への写像を$f:\mathcal Z \rightarrow \mathcal X=\mathbb{R}^D$とし，$\mathbf x\in\mathcal X$は次元に関して独立と仮定します．　すなわち，

\begin{align}
x_{nd}=f_d(\mathbf z_n)+\epsilon_{nd}
\end{align}

にしたがって$\mathbf z_n$から$x_{nd}$がそれぞれの次元ごとに独立に生成されると仮定します．ここで$\epsilon_{nd}$は精度パラメータを$\beta$とするガウスノイズです．$f_d$の事前分布を$f_d\sim \mathcal{GP}(0,k(\mathbf{z},\mathbf{z}'))$とします¹¹．観測データを$\mathbf X =(\mathbf x_1,\mathbf x_2,\cdots,\mathbf x_N)$とし対応する潜在変数を$\mathbf Z$とした時，事前分布は$p(\mathbf f_d\mid \mathbf Z)=\mathcal{N}(\mathbf f_d\mid \mathbf 0, k(\mathbf Z,\mathbf Z))$となります¹²．ここで$\mathbf f_d=(f_d(\mathbf z_1),f_d(\mathbf z_2),\cdots,f_d(\mathbf z_N))$です．この時$\mathbf f_d$は次元$d$について独立なので周辺尤度は

\begin{align}
p(\mathbf{X}\mid \mathbf{Z},\beta) &= \prod^D_{d=1}p(\mathbf{x}_{:d}\mid \mathbf{Z},\beta) \\
&= \prod^D_{d=1}\int{p(\mathbf{x}_{:d}\mid \mathbf{f}_d,\beta)p(\mathbf{f}_d\mid \mathbf{Z})}d\mathbf{f}_d
\end{align}

となります．また$p(\mathbf{x}_{:d}\mid \mathbf{f}_d,\beta)$も$p(\mathbf{f}_d\mid \mathbf{Z})$もガウス分布となるので，周辺尤度もガウス分布となります．

したがって$p(\mathbf x_{:d}\mid\mathbf Z)=\mathcal N(\mathbf x_{:d}\mid\mathbb E[\mathbf x_{:d}],\mathbb V[\mathbf x_{:d}])$とすると以下のようになります．

\begin{align}
\mathbb{E}[\mathbf x_{:d}]&=\mathbb{E}[\mathbf f_d+\boldsymbol{\epsilon}_{:d}]\\
&=\mathbf 0\\
\mathbb{V}[\mathbf x_{:d}]&=\mathbb{E}[(\mathbf f_d+\boldsymbol{\epsilon}_{:d})(\mathbf f_d+\boldsymbol{\epsilon}_{:d})^{\rm T}] \\
&=\mathbb{E}[\mathbf f_d\mathbf f^{\rm T}_d]+2\mathbb{E}[\boldsymbol{\epsilon}_{:d}\mathbf f^{\rm T}_d]+\mathbb{E}[\boldsymbol{\epsilon}_{:d}\boldsymbol{\epsilon}^{\rm T}_{:d}] \\
&=k(\mathbf Z,\mathbf Z)+\beta^{-1}\mathbf I_N \\
&=\mathbf K
\end{align}

これより

\begin{align}
\log{p(\mathbf{X}\mid \mathbf{Z},\beta)} &= \sum^D_{d=1}\log{p(\mathbf{x}_{:d}\mid \mathbf{Z},\beta)} \\
&= -\frac{ND}{2}\log{2\pi}-\frac{D}{2}\log{|\mathbf{K}|}-\frac{D}{2}{\rm Tr}[\mathbf{K}^{-1}\mathbf{S}]\\
\end{align}

となり，式(5)と一致します．　つまり， DPPCAで推定する潜在変数$\mathbf Z$は見方を変えると，入力が潜在変数になる代わりに出力が多次元のガウス過程を考えた時に周辺尤度を最大にする潜在変数を推定することと等価であることがわかります．このことから潜在空間から観測空間への写像$f$の事後分布は以下のようなガウス過程として書くことができます．

\begin{align}
f_d(\mathbf{z}) &\sim \mathcal{GP}(\mu_d(\mathbf{z}),\sigma(\mathbf{z},\mathbf{z}')) \\
\mu_d(\mathbf{z}) &= k(\mathbf{z},\mathbf{Z}) (k(\mathbf{Z},\mathbf{Z})+\beta^{-1}\mathbf{I})^{-1}\mathbf{x}_{:d} \\
\sigma(\mathbf{z},\mathbf{z}') &= k(\mathbf{z},\mathbf{z}') - k(\mathbf{z},\mathbf{Z})(k(\mathbf{Z},\mathbf{Z})+\beta^{-1}\mathbf{I})^{-1}k(\mathbf{Z},\mathbf{z}')
\end{align}

さらにハイパーパラメータについても式(5)を最大にするハイパーパラメータを推定することが周辺尤度を最大にするハイパーパラメータを推定することと等価であることがわかるため教師なし学習におけるハイパーパラメータ推定の妥当性も周辺尤度を最大化しているという観点で保証することができます．

式(5)は解析的に解くことができないので勾配法を用いて$\mathbf Z$とハイパーパラメータの推定を行います．$z_{nl}$について$L_{\rm DPPCA}$を微分する際，$z_{nl}$はカーネル関数の引数となっているため，以下の微分の連鎖律を用いて計算することができます．

\begin{align}
\frac{\partial L_{\rm DPPCA}}{\partial z_{nl}}={\rm Tr}\left(\frac{\partial L_{\rm DPPCA}}{\partial \mathbf{K}}\frac{\partial \mathbf{K}}{\partial z_{nl}}\right)
\end{align}

ここで

\begin{align}
\frac{\partial L_{\rm DPPCA}}{\partial \mathbf{K}} = \frac{D}{2}(\mathbf{K}^{-1}\mathbf{S}\mathbf{K}^{-1}-\mathbf{K}^{-1})
\end{align}

となります．$\frac{\partial \mathbf{K}}{\partial z_{nl}}$についてはカーネル関数に依存するためそれに合わせて微分してください．また，実際のアルゴリズムの計算では潜在変数が極端な値になることを抑えるために事前分布の対数$\log{p(\mathbf{Z})}$を加えた式を最大化します．事前分布は基本は標準正規分布が用いられます．

GPLVMの動作検証

最後に実際にGPLVMをpython上で実装して簡単なデータで学習がうまく行くかどうか検証します．
今回は潜在変数$\mathbf z$を$[-1,1]^2$の範囲で一様に200点ランダムサンプリングしたデータを以下の関数で３次元空間へ写像したデータを観測データとしました．ここで$\mathbf{\epsilon}$はガウスノイズです．またカーネル関数はガウスカーネルを使い，ハイパーパラメータと観測ノイズも周辺尤度最大化で推定しています．潜在変数の初期値もランダムに決めています．

\begin{align}
x_{n1} &= z_{n1}+\epsilon_{n1} \\
x_{n2} &= z_{n2}+\epsilon_{n2} \\
x_{n3} &= z_{n1}^2 - z_{n2}^2+\epsilon_{n3} \\
\end{align}

実際のプログラムは以下のようになっています．

GPLVM.py

import numpy as np

class GPLVM(object):
    def __init__(self,Y,LatentDim,HyperParam,X=None):
        self.Y = Y
        self.hyperparam = HyperParam
        self.dataNum = self.Y.shape[0]
        self.dataDim = self.Y.shape[1]

        self.latentDim = LatentDim
        if X is not None:
            self.X = X
        else:
            self.X = 0.1*np.random.randn(self.dataNum,self.latentDim)
        self.S = Y @ Y.T
        self.history = {}

    def fit(self,epoch=100,epsilonX=0.5,epsilonSigma=0.0025,epsilonAlpha=0.00005):

        resolution = 10
        M = resolution**self.latentDim
        self.history['X'] = np.zeros((epoch, self.dataNum, self.latentDim))
        self.history['F'] = np.zeros((epoch, M, self.dataDim))
        sigma = np.log(self.hyperparam[0])
        alpha = np.log(self.hyperparam[1])
        for i in range(epoch):

            # 潜在変数の更新
            K = self.kernel(self.X,self.X,self.hyperparam[0]) + self.hyperparam[1]*np.eye(self.dataNum)
            Kinv = np.linalg.inv(K)
            G = 0.5*(Kinv @ self.S @ Kinv-self.dataDim*Kinv)
            dKdX = -2.0*(((self.X[:,None,:]-self.X[None,:,:])*K[:,:,None]))/self.hyperparam[0]
            # dFdX = (G[:,:,None] * dKdX).sum(axis=1)-self.X
            dFdX = (G[:,:,None] * dKdX).sum(axis=1)

            # ハイパーパラメータの更新
            Dist = ((self.X[:, None, :] - self.X[None, :, :]) ** 2).sum(axis=2)
            dKdSigma = 0.5*Dist/self.hyperparam[0] * ( K- self.hyperparam[1]*np.eye(self.dataNum) )
            dFdSigma = np.trace(G @ dKdSigma)

            dKdAlpha = self.hyperparam[1]*np.eye(self.dataNum)
            dFdAlpha = np.trace(G @ dKdAlpha)

            self.X = self.X + epsilonX * dFdX
            self.history['X'][i] = self.X
            sigma = sigma + epsilonSigma * dFdSigma
            self.hyperparam[0] = np.exp(sigma)
            alpha = alpha + epsilonAlpha * dFdAlpha
            self.hyperparam[1] = np.exp(alpha)

            zeta = np.meshgrid(np.linspace(self.X[:, 0].min(), self.X[:, 0].max(), resolution),
                               np.linspace(self.X[:, 1].min(), self.X[:, 1].max(), resolution))
            zeta = np.dstack(zeta).reshape(M, self.latentDim)
            K = self.kernel(self.X,self.X,self.hyperparam[0]) + self.hyperparam[1]*np.eye(self.dataNum)
            Kinv = np.linalg.inv(K)
            KStar = self.kernel(zeta, self.X, self.hyperparam[0])
            self.F = KStar @ Kinv @ self.Y
            self.history['F'][i] = self.F


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

main.py

from GPLVM import GPLVM
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import os

def createKuraData(N, D,sigma=0.1):
    X = (np.random.rand(N, 2) - 0.5) * 2
    Y = np.zeros((N, D))
    Y[:, :2] = X
    Y[:,2]=X[:,0]**2-X[:,1]**2
    Y += np.random.normal(0,sigma,(N,D))

    return [X,Y]

def plot_prediction(Y,f,Z,epoch,save_folder):
    fig = plt.figure(1,[10,8])

    nb_nodes = f.shape[1]
    nb_dim = f.shape[2]
    resolution = np.sqrt(nb_nodes).astype('int')
    for i in range(epoch):
        if i%10 is 0:
            ax_input = fig.add_subplot(1, 2, 1, projection='3d')
            ax_input.cla()
            ax_input.set_aspect("auto")
            # 観測空間の表示
            r_f = f[i].reshape(resolution, resolution, nb_dim)
            ax_input.plot_wireframe(r_f[:, :, 0], r_f[:, :, 1], r_f[:, :, 2],color='k')
            ax_input.scatter(Y[:, 0], Y[:, 1], Y[:, 2], c=Y[:, 0], edgecolors="k",marker='x')
            ax_input.set_xlim(Y[:, 0].min(), Y[:, 0].max())
            ax_input.set_ylim(Y[:, 1].min(), Y[:, 1].max())
            ax_input.set_zlim(Y[:, 2].min(), Y[:, 2].max())

            # 潜在空間の表示
            ax_latent = fig.add_subplot(1, 2, 2)
            ax_latent.cla()
            ax_latent.set_aspect("equal")
            ax_latent.set_xlim(Z[:,:, 0].min(), Z[:,:, 0].max())
            ax_latent.set_ylim(Z[:,:, 1].min(), Z[:,:, 1].max())
            ax_latent.scatter(Z[i,:, 0], Z[i,:, 1], c=Y[:, 0], edgecolors="k")
            plt.savefig(save_folder+"fig{0}.png".format(i))

        plt.pause(0.001)
    plt.show()

if __name__ == '__main__':
    L=2
    N=200
    D=3
    sigma=3
    alpha=0.05
    beta=0.08
    seedData=1
    resolution = 10
    M = resolution**L
    save_folder = "fig/"
    os.makedirs(save_folder,exist_ok=True)

    # 入力データの生成
    # np.random.seed(seedData)
    [X,Y] = createKuraData(N,D,sigma=0.01)

    # カーネルの設定
    [U,D,Vt] = np.linalg.svd(Y)
    model = GPLVM(Y,L, np.array([sigma**2,alpha/beta]))
    # GPLVMの最適化
    epoch = 200
    model.fit(epoch=epoch,epsilonX=0.05,epsilonSigma=0.0005,epsilonAlpha=0.00001)

    # 推定した潜在変数の取得
    X = model.history['X']
    f = model.history['F']

    # 学習結果の表示
    plot_prediction(Y,f,X,epoch,save_folder)

左図が観測空間上での推定した多様体の学習過程です．右図はその時の潜在変数です．パラメータの初期値がある程度悪くなければ初期値をランダムに決めても安定して動作しています．特にカーネル関数の幅をあらかじめ大きめにとっておくと安定しやすい傾向があります．　

おわりに

私のGPLVMに対する理解をまとめてみました．おかしい部分やご助言等ございましたらご連絡いただけると幸いです．

Jack Wang and Hertzmann, Aaron and Fleet, David J,　Gaussian Process Dynamical Models，Advances in Neural Information Processing Systems 18,2006. ↩
Carl Henrik. Ek, Shared Gaussian Process Latent Variable Models, PhD Thesis, 2009. ↩
Zenglin Xu and Feng Yan and Yuan Qi, Infinite Tucker Decomposition: Nonparametric Bayesian Models for Multiway Data Analysis, Proceedings of the 29th International Conference on Machine Learning, ICML, 2012. ↩
Neil Lawrence, Probabilistic Non-linear Principal Component Analysis with Gaussian Process Latent Variable Models, J. Mach. Learn. Res., 2005. ↩
今回は$\mathbf{x}$の平均がゼロであると仮定しますが本質的には違いはありません． ↩
パッと見た感じではPCAとpPCAの目的関数は異なっているように見え，PCAとpPCAの目的関数はどれくらい類似しているのか疑問です． ↩
Michael Tipping and Christopher Bishop, Probabilistic Principal Component Analysis, Journal of the Royal Statistical Society Series, 1999. ↩
厳密にはカーネル主成分分析はDPPCAではなく通常の主成分分析の双対表現を用いるため異なりますがここではそこの違いは無視しています． ↩
James T. Kwok and Ivor W. Tsang, The pre-image problem in kernel methods, , IEEE Transactions on Neural Networks, 2004. ↩
Carl E. Rasmussen and Christopher K. I. Williams. Gaussian Processes for Machine Learning. The MIT Press, 2006. ↩
$\mathcal{GP}(\mu_0(\mathbf{z}),k(\mathbf{z},\mathbf{z}))$は平均関数$\mu_0(\mathbf{z})$，共分散関数$k(\mathbf{z},\mathbf{z}')$のガウス過程です． ↩
ガウス過程に任意の入力集合を与えたときの出力の確率分布は無矛盾にガウス分布に従う性質を利用しています． ↩

GPLVMの基本形についてまとめてみた