ガウス過程力学モデル(Gaussian Process Dynamical Model,以下GPDM) はガウス過程による教師なし学習の一つで、 ガウス過程潜在モデル(Gaussian Process Latent Variable Model,以下GPLVM) を拡張したモデルのうちの一つです。

GPLVMでは、潜在変数 $\mathbf{X}=(\mathbf{x}_1,\mathbf{x}_2,\dots,\mathbf{x}_N)$ の独立性を仮定していました。これに対してGPDMは、潜在変数 $\mathbf{X}$ が時系列データであるという仮定を導入して、潜在空間での構造を学習します。

この記事では、GPDMの原論文を参考にしながら、GPDMをゼロから実装し、PCAやGPLVMなどの他の次元圧縮手法と比較してみます。GPDMの基礎となるGPLVMに関する解説は行いません。GPLVMについても知りたい方はこちらの記事が参考になると思います。

※注意
この記事は筆者が勉強したことをまとめて、理解を深めるために執筆しています。内容に誤り等があった場合は、コメント等でご指摘いただけると幸いです。

ガウス過程力学モデル(Gaussian Process Dynamical Model)

このセクションでは、時系列性のある観測 $\mathbf{Y}$ が得られたときの潜在変数 $\mathbf{X}$ の分布 $p(\mathbf{X}|\mathbf{Y})$ をベイズの定理を用いて導出していきます。

p(\mathbf{X}|\mathbf{Y})\propto p(\mathbf{Y}|\mathbf{X})p(\mathbf{X})

J次元の観測データ $\mathbf{Y}=(\mathbf{y}_1,\mathbf{y}_2,\dots,\mathbf{y}_T)$ とD次元の潜在変数 $\mathbf{X}=(\mathbf{x}_1,\mathbf{x}_2,\dots,\mathbf{x}_T)$ に対してマルコフ性を仮定します。

\mathbf{x}_t=f(\mathbf{x}_{t-1};\mathbf{A})+\mathbf{n}_{x,t} \quad \mathbf{x}_t\in\mathbb{R}^J \quad

\mathbf{y}_t=g(\mathbf{x}_{t};\mathbf{B})+\mathbf{n}_{y,t} \quad \mathbf{y}_t\in\mathbb{R}^D

確率モデルとして $\mathbf{f}$ と $\mathbf{g}$ に対して平均ゼロのガウシアンノイズが加わったものを考えます。 $f$ と $g$ は次の形で与えられます。

f(\mathbf{x};\mathbf{A})=\sum_{k=1}^{K}\mathbf{a}_{k} \phi_k(\mathbf{x})=(\mathbf{a}_1,\mathbf{a}_2,\dots,\mathbf{a}_K)
\left(
\begin{array}{c}
\phi_1(\mathbf{x}) \\
\phi_2(\mathbf{x}) \\
\vdots \\
\phi_K(\mathbf{x})
\end{array}
\right)

g(\mathbf{x};\mathbf{B})=\sum_{m=1}^{M}\mathbf{b}_{m} \psi_m(\mathbf{x})=(\mathbf{b}_1,\mathbf{b}_2,\dots,\mathbf{b}_M)
\left(
\begin{array}{c}
\psi_1(\mathbf{x}) \\
\psi_2(\mathbf{x}) \\
\vdots \\
\psi_M(\mathbf{x})
\end{array}
\right)

原論文では、 $\mathbf{B}$ の各列ベクトル $\mathbf{b}_j$ が独立にガウス分布に従うと仮定しています。すなわち、分散を $w_j^{-2}$ とすると

p(\mathbf{B})=\prod_{j=1}^{J}\mathbb{N}(\mathbf{b}_j|\mathbf{0},w_j^{-2}\mathbf{I})

と記述することができます。

論文の仮定に従い $\mathbf{n}_{y,j}$ の分散を $w_j^{-2}\sigma_Y^2$ とすると、$\mathbf{b}_j$ が与えられたときの $\mathbf{y}_j$ は $T$ 次元のガウス分布に従います。

\mathbf{y}_j=\mathbf{\Psi}\mathbf{b}_j+\mathbf{n}_{y,j}

\mathbf{\Psi}=(\psi(\mathbf{x}_1),\psi(\mathbf{x}_2),\dots,\psi(\mathbf{x}_T))^{\top}

\mathbf{y}_j|\mathbf{X},\mathbf{b}_j\sim\mathbb{N}_T(\mathbf{y}_j|\mathbf{\Psi}\mathbf{b}_j,w_j^{-2}\sigma_Y^2\mathbf{I})

さらに $p(\mathbf{b}_j)$ と $p(\mathbf{y}_j|\mathbf{b}_j)$ が共に独立にガウス分布に従うことから、 $p(\mathbf{y}_j|\mathbf{b}_j)$ を　$\mathbf{b}_j$ について周辺化することができます。

\begin{align}
p(\mathbf{y}_j|\mathbf{X})&=\mathbb{N}_T(\mathbf{y}_j|\mathbf{0},\mathbf{\Psi}w_j^{-2}\mathbf{\Psi}^{\top}+w_j^{-2}\sigma_Y^{2}\mathbf{I})\\
&=\mathbb{N}_T(\mathbf{y}_j|\mathbf{0},w_j^{-2}(\mathbf{\Psi}\mathbf{\Psi}^{\top}+\sigma_Y^{2}\mathbf{I}))\\
&=\mathbb{N}_T(\mathbf{y}_j|\mathbf{0},w_j^{-2}\mathbf{K}_Y)
\end{align}

\mathbf{K}_Y=w_j^{-2}(\mathbf{\Psi}\mathbf{\Psi}^{\top}+\sigma_Y^{2}\mathbf{I})

潜在変数 $\mathbf{X}$ から各観測 $\mathbf{y}_j$ が生成される確率 $p(\mathbf{y}_j|\mathbf{X})$ を求められたので、これを用いて $p(\mathbf{Y}|\mathbf{X})$ を求めると

\begin{align}
p(\mathbf{Y}|\mathbf{X})&=\prod_{j=1}^{j}p(\mathbf{y}_j|\mathbf{X})\\
&=\prod_{j=1}^{j}\frac{1}{\sqrt{(2\pi)^T |w_j^{-2}\mathbf{K}_Y|}}\exp\{-\frac{1}{2}\mathbf{y}_j^{\top}w_j^{-2}\mathbf{K}_Y^{-1}\mathbf{y}_j\}\\
&=\frac{(w_j^{-2})^{TJ}}{\sqrt{(2\pi)^{TJ} |\mathbf{K}_Y|^{J}}}\exp\{-\frac{1}{2}\sum_{j=1}^{J}\mathbf{y}_j^{\top}w_j^{-2}\mathbf{K}_Y^{-1}\mathbf{y}_j\}\\
&=\frac{(w_j^{-2})^{TJ}}{\sqrt{(2\pi)^{TJ} |\mathbf{K}_Y|^{J}}}\exp\{-\frac{1}{2}\textrm{tr}(\sum_{j=1}^{J}\mathbf{K}_Y^{-1}\mathbf{y}_j\mathbf{y}^{\top}_jw_j^{-2})\}\\
&=\frac{|\mathbf{W}|^{T}}{\sqrt{(2\pi)^{TJ} |\mathbf{K}_Y|^{J}}}\exp\{-\frac{1}{2}\textrm{tr}(\mathbf{K}_Y^{-1}\mathbf{Y}\mathbf{W}\mathbf{Y}^{\top})\}
\end{align}

\mathbf{W}=\textrm{diag}([w_1^{-2},w_2^{-2},\dots,w_J^{-2}])

が得られます。

同様に $p(\mathbf{X})$ を $\mathbf{A}$ の各列ベクトル $\mathbf{a}_d$ が独立にガウス分布に従うと仮定して求めていきます。

p(\mathbf{A})=\prod_{j=d}^{D}\mathbb{N}_K(\mathbf{a}_d|\mathbf{0},\mathbf{I})

\mathbf{x}_d^{(2:T)}=\mathbf{\Phi}^{1:(T-1)}\mathbf{a}_d+\mathbf{n}_{x,d}

\mathbf{x}_d^{(2:T)}|\mathbf{X},\mathbf{a}_d\sim\mathbb{N}_{T-1}(\mathbf{x}_d^{(2:T)}|\mathbf{\Phi}^{1:(T-1)}\mathbf{a}_d,\sigma_X^2\mathbf{I})

$\mathbf{x}_d^{(2:T)}$ は $\mathbf{X}$ の第1行を $\mathbf{X}$ から取り除いた　$\mathbf{X}$ のd番目の列のベクトル $\mathbf{x}_d^{(2:T)}=(x_d^{(2)},\dots,x_d^{(T)})^\top$ を表しています。また、 $\mathbf{\Phi}^{1:(T-1)}$ は第T行を　$\mathbf{\Phi}$　から取り除いたベクトル $\mathbf{\Phi}^{1:(T-1)}=(\phi(\mathbf{x}_1),\phi(\mathbf{x}_2),\dots,\phi(\mathbf{x} _{T-1}))^{\top}$ に対応しています。

$p(\mathbf{x}_d^{(2:T)}|\mathbf{a}_d)$ が得られたので、周辺化をおこなって $p(\mathbf{x}_d)$ を求めます。

\begin{align}
p(\mathbf{x}_d)&=\mathbb{N}_{T-1}(\mathbf{x}_d|\mathbf{0},\mathbf{\Phi}\mathbf{\Phi}^{\top}+\sigma_X^{2}\mathbf{I})\\
&=\mathbb{N}_{T-1}(\mathbf{x}_d|\mathbf{0},\mathbf{K}_X)
\end{align}

\mathbf{K}_X=\mathbf{\Phi}\mathbf{\Phi}^{\top}+\sigma_X^{2}\mathbf{I}

以上の手順より、 $\mathbf{X}$ の分布は次のように求められます。

\begin{align}
p(\mathbf{X})&=\int p(\mathbf{X}|\mathbf{A})p(\mathbf{A})d\mathbf{A}\\
&=\int p(\mathbf{x}_1)\prod_{d=1}^{D}p(\mathbf{x}_d|\mathbf{A})p(\mathbf{A})d\mathbf{A}\\
&=p(\mathbf{x}_1)\prod_{d=1}^{D}\frac{1}{\sqrt{(2\pi)^{T-1} |\mathbf{K}_X|}}\exp\{-\frac{1}{2}\mathbf{x}_d^{\top}\mathbf{K}_X^{-1}\mathbf{x}_d\}\\
&=p(\mathbf{x}_1)\frac{1}{\sqrt{(2\pi)^{D(T-1)} |\mathbf{K}_X|^{D}}}\exp\{-\frac{1}{2}\sum_{d=1}^{D}\mathbf{x}_d^{\top}\mathbf{K}_X^{-1}\mathbf{x}_d\}\\
&=p(\mathbf{x}_1)\frac{1}{\sqrt{(2\pi)^{D(T-1)} |\mathbf{K}_X|^{D}}}\exp\{-\frac{1}{2}\textrm{tr}(\mathbf{K}_X^{-1}\mathbf{X}^{(2:T)}\mathbf{X}^{(2:T)\top})\}
\end{align}

時系列性のある観測 $\mathbf{Y}$ が得られたときの潜在変数 $\mathbf{X}$ の分布 $p(\mathbf{X}|\mathbf{Y})$ はベイズの定理 $p(\mathbf{X}|\mathbf{Y})\propto p(\mathbf{Y}|\mathbf{X})p(\mathbf{X})$ を用いると

\begin{align}
\log  p(\mathbf{X}|\mathbf{Y})&=\log p(\mathbf{Y}|\mathbf{X})+p(\mathbf{X})\\
&=T\log |\mathbf{W}|-\frac{TJ}{2}\log (2\pi)-\frac{J}{2}\log |\mathbf{K}_X|-\frac{1}{2}\textrm{tr}
(\mathbf{K}_Y^{-1}\mathbf{Y}\mathbf{W}\mathbf{Y}^{\top})\\
&+\frac{D(T-1)}{2}\log(2\pi)-\frac{D}{2}\log |\mathbf{K}_X|-\frac{1}{2}\textrm{tr}(\mathbf{K}_X^{-1}\mathbf{X}^{(2:T)}\mathbf{X}^{(2:T)\top})+\log p(\mathbf{x}_1)\\
&=T\log |\mathbf{W}|-\frac{J}{2}\log |\mathbf{K}_X|-\frac{1}{2}\textrm{tr}
(\mathbf{K}_Y^{-1}\mathbf{Y}\mathbf{W}\mathbf{Y}^{\top})-\frac{D}{2}\log |\mathbf{K}_X|-\frac{1}{2}\textrm{tr}(\mathbf{K}_X^{-1}\mathbf{X}^{(2:T)}\mathbf{X}^{(2:T)\top})
\end{align}

となります。上記の計算では簡単化のために対数尤度を計算しました。

なお、二つのカーネル行列 $\mathbf{K}_X$ と $\mathbf{K}_Y$ に用いられるカーネルは、論文では次のように提示されています。先ほど分散として用いた $\sigma_X^2$ と　$\sigma_Y^2$ はそれぞれ $\beta_3^{-1}$ と $\alpha_4^{-1}$ に対応していることがわかります。

k_Y(\mathbf{x},\mathbf{x}')=\beta_1\exp(-\frac{\beta_2}{2}\|\mathbf{x}-\mathbf{x}'\|)+\beta_3^{-1}\delta_{\mathbf{x},\mathbf{x}'}

k_X(\mathbf{x},\mathbf{x}')=\alpha_1\exp(-\frac{\alpha_2}{2}\|\mathbf{x}-\mathbf{x}'\|)+\alpha_3\mathbf{x}^{\top}\mathbf{x}'+\alpha_4^{-1}\delta_{\mathbf{x},\mathbf{x}'}

実装

GPDMをsklearnのS字カーブのデータに対して用いてみます。sklearnのmake_s_curve()は元々3次元のデータセットを返しますが、扱いやすくするために軸を一つ削除して2次元のデータに変換します。

S字カーブのデータからカーネル行列を定め、平均0で正規分布に従った40次元のデータにノイズをのせたデータセットを作成します。

この実験では、40次元の観測データ $\mathbf{Y}$ から2次元の潜在変数 $\mathbf{X}$ をGPDMで推定し、PCAやGPLVMなどの他の次元圧縮手法と比較検討を行います。

import autograd.numpy as np
from sklearn.datasets import make_s_curve

r = np.random.RandomState(42)

def make_dataset(n_samples=200, n_dims=40):
    X, t = make_s_curve(n_samples, random_state=42)

    X = np.delete(X, obj=1, axis=1)

    indices = t.argsort()
    X, t = X[indices], t[indices]

    K = kernel(X, 1, 2)

    F = r.multivariate_normal(np.zeros(n_samples), K, size=n_dims).T
    Y = F + r.normal(0, scale=1, size=F.shape)
    return X, Y, t


def rbf(x, x_prime, theta_1, theta_2):

    return theta_1 * np.exp(-1 * np.sum((x - x_prime)**2) / theta_2)


def kernel(X, theta_1, theta_2):
    length = X.shape[0]

    K = np.zeros((length, length))
    for x_idx in range(length):
        for x_prime_idx in range(length):
            K[x_idx, x_prime_idx] += rbf(X[x_idx], X[x_prime_idx], theta_1, theta_2)

    return K

次に対数尤度を計算するプログラムを実装します。基本的には一つ前のセクションで導出した式をそのままプログラムに落とし込むだけです。

def kernel_Y(X, theta_1, theta_2, theta_3):
    length = X.shape[0]
    diffs = np.expand_dims(X / theta_2, 1) - np.expand_dims(X / theta_2, 0)
    return theta_1 * np.exp(-0.5 * np.sum(diffs ** 2, axis=2)) + theta_3 * np.eye(length)


def kernel_X(X, theta_1, theta_2, theta_3, theta_4):
    K = kernel_Y(X, theta_1, theta_2, theta_3)
    return K + theta_4 * X @ X.T


def log_posterior(Y, X, beta, alpha):
    _, n_dims = Y.shape
    D = X.shape[1]

    K_Y = kernel_Y(X, *beta)
    det_term = -n_dims * np.prod(np.linalg.slogdet(K_Y)) / 2
    tr_term = -1 * np.trace(np.linalg.inv(K_Y) @ Y @ Y.T) / 2
    LY = det_term + tr_term

    K_X = kernel_X(X[:-1], *alpha)
    x = X[1:]
    det_term = - D * np.prod(np.linalg.slogdet(K_X)) / 2
    tr_term = -1 * np.trace(np.linalg.inv(K_X) @ x @ x.T) / 2
    LX = det_term + tr_term
    return - LY - LX

kernel_Yの実装が少しわかりずらいものなっていますが、上の実装は下の実装と等価になっています。なぜ、np.expand_dims()などの関数を用いて記述したのかというと、勾配の計算を省略するためにautogradでnumpyをラップしているためです。

解析的に関数の勾配を求める作業は手間がかかることが多く、勾配の式をそのままプログラムに落とし込む際にも勾配の式が長ければ長いほど実装ミスが起こる可能性も増えてしまいます。autogradでラップしてnumpyをインポートすると通常のnumpyの関数が使える上にgrad()というautogradの関数に関数を渡すことで自動で勾配を計算してくれます。

しかしながら、勾配計算中はArrayBoxと呼ばれるautogradのクラスに値が格納されるため、下のような実装をしてしまうとX[x_idx]などで値にアクセスできません。(X._values[x_idx]で値にアクセスできます。) したがって、autogradを用いるときはなるべくインデックスによるアクセスは避けて計算を行う必要があります。

def rbf(x, x_prime, theta_1, theta_2):

    return theta_1 * np.exp(-1 * np.sum((x - x_prime)**2) / theta_2)


# Radiant Basis Kernel + Error
def kernel_Y(X, theta_1, theta_2, theta_3):
    length = X.shape[0]

    K = np.zeros((length, length))
    for x_idx in range(length):
        for x_prime_idx in range(length):
            k = rbf(X[x_idx], X[x_prime_idx], theta_1, theta_2)
            K[x_idx, x_prime_idx] += k

    K += theta_3 * np.eye(length)
    return K

最後に勾配法によってハイパーパラメータや潜在変数 $\mathbf{X}$ を最適化するプログラムを実装します。初期値としては、完全にランダムな初期値でも問題ないですが、GPyなどのライブラリでGPLVMがPCAで初期化をおこなっていたことに倣い、PCAで初期化をおこなっています。

勾配法のアルゴリズムにはL-BFGS-B法を用いました。こちらも勾配法であれば性能差こそあれ、どんなアルゴリズムでも問題ありません。

from autograd import grad
from scipy.optimize import fmin_l_bfgs_b
from sklearn.decomposition import PCA

def optimize_gpdm(Y, n_components):
    n_samples = Y.shape[0]

    pca = PCA(n_components=n_components)
    X_init = pca.fit_transform(Y)
    beta0 = np.array([1, 1, 1e-6])
    alpha0 = np.array([1, 1, 1e-6, 1e-6])

    def _lml(params):
        X = params[:n_samples * n_components].reshape(X_init.shape)
        beta = params[n_samples * n_components:n_samples * n_components + 3]
        alpha = params[n_samples * n_components + 3:]
        return log_posterior(Y, X, beta, alpha)

    def obj_func(params):
        _grad = grad(_lml)
        return _lml(params), _grad(params)

    x0 = np.concatenate([X_init.flatten(), beta0, alpha0])
    opt_res = fmin_l_bfgs_b(obj_func, x0, epsilon=1e-32)
    X_map = opt_res[0][:n_samples * n_components].reshape(X_init.shape)
    return X_map

GPDMを用いて次元圧縮するプログラムを実行してみます。比較のためPCAとGPLVMによる次元圧縮も示しました。

import matplotlib.pyplot as plt
from sklearn.decomposition import PCA
from GPy.models import GPLVM
plt.style.use('seaborn-pastel')

def plot(X_map, title):
    fig = plt.figure(figsize=(16, 8))
    fig.suptitle(f"{title}", fontsize=18, fontweight='bold')

    ax1 = fig.add_subplot(1, 2, 1)
    ax1.set_title("Scatter", fontsize=18)
    ax1.scatter(X_map[:, 0], X_map[:, 1], c=t, cmap='Blues')

    ax2 = fig.add_subplot(1, 2, 2)
    ax2.set_title("Line", fontsize=18)
    ax2.plot(X_map[:, 0], X_map[:, 1])

    fig.tight_layout()
    fig.savefig("{}.png".format(title.lower()))
    plt.close()


if __name__ == "__main__":
    X, Y, t = make_dataset(n_samples=200)

    # original data
    fig = plt.figure()
    plt.scatter(X[:, 0], X[:, 1], c=t, cmap='Blues')
    fig.savefig("scurve.png")

    # pca
    pca = PCA(n_components=2)
    X_map = pca.fit_transform(Y)
    plot(X_map, title='PCA')

    # gplvm
    gplvm = GPLVM(Y, input_dim=2)
    gplvm.optimize()
    X_map = gplvm.X
    plot(X_map, title='GPLVM')

    # gpdm
    X_map = optimize_gpdm(Y, n_components=2)
    plot(X_map, title='GPDM')

一つ目の結果はPCAによる次元圧縮の結果です。40次元のデータを2次元に圧縮しています。潜在変数であるS字カーブがそもそも多峰性のデータであるので、PCAでは潜在変数をうまく推定できていません。

二つ目の結果はGPLVMによる次元圧縮の結果です。PCAよりもラインプロットが多少滑らかになっているように見えます。また、散布図もS字を描いているようにも見えなくもないです。

三つ目の結果はGPDMによる次元圧縮の結果です。GPLVMとは違いデータに時系列性を過程しているので、今回のデータセットに対する次元圧縮手法としては3つの中で理論的に最も適しています。実際に散布図もラインプロットも3つの次元圧縮アルゴリズムの中で最も滑らかかつS字のカーブを推定できています。

まとめ

GPDMをゼロから実装し、PCAやGPDMなどの他の次元圧縮アルゴリズムと比較してみました。データに時系列性があることがわかっているときは、次元圧縮手法として有用なアルゴリズムと言えそうです。

今回の実験コードはGitHub上にあります。

参考

[1] 持橋大地, ガウス過程と機械学習. 2019.
[2] Jack M. Wang, David J. Fleet, Aaron Hertzmann, Gaussian Process Dynamical Models. 2008.
[3] Gregory Gundersen, Gaussian Process Dynamical Models. 2020.