[PYTHON] J'ai essayé de résumer la forme de base de GPLVM

Cet article est tiré du Furukawa Lab Advent_calendar Day 23.

introduction

GPLVM est une méthode d'apprentissage non supervisée qui utilise le processus gaussien (GP) et peut estimer des variantes de faible dimension dans lesquelles les données sont distribuées dans un espace de grande dimension. Le point le plus intéressant de la GPLVM est que le modèle est très simple, il est donc très extensible et théoriquement facile à manipuler. En fait, il a été étendu à des situations plus complexes telles que l'estimation d'hyperparamètres dans le cadre de l'estimation bayésienne, l'analyse des données de séries chronologiques [^ 1], l'analyse multivue [^ 2] et la décomposition tenseur [^ 3]. .. Cette fois, l'algorithme GPLVM le plus basique est dérivé en partant de l'ACP probabiliste (pPCA). Cette dérivation suit généralement le flux de l'article du professeur Lawrence, veuillez donc vous référer à cette personne pour plus de détails [^ 4].

probabilistic PCA

Comme son nom l'indique, le pPCA est un modèle qui redistribue l'analyse en composantes principales dans le cadre de la théorie des probabilités. Dans pPCA, $ \ mathbf {X} = (\ mathbf {x} _1, \ mathbf {x} _2, \ cdots, \ mathbf {x} _N) \ in \ mathbb {R} ^ {D \ times N} $ Supposons que $ \ mathbf {x} $ soit observé en mappant la variable latente $ \ mathbf {z} \ in \ mathbb {R} ^ L $ donnée comme suit [^ 5] ].

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

Ici, $ \ epsilon $ est le bruit observé, et le bruit gaussien isotrope avec le paramètre de précision $ \ beta ^ {-1} $ est supposé. C'est aussi $ \ mathbf {W} \ in \ mathbb {R} ^ {D \ times L} $. Nous supposons également que $ \ mathbf {z} $ suit une distribution normale standard de la dimension $ L $. Autrement dit, $ \ mathbf {x} $ et $ \ mathbf {z} $ suivent la distribution de probabilité suivante.

\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}

À ce stade, $ p (\ mathbf {x} \ mid \ mathbf {W}, \ mathbf {z}, \ beta) $ et $ p (\ mathbf {z}) $ sont des distributions gaussiennes, donc distribution simultanée $ p (\ mathbf {x}, \ mathbf {z} \ mid \ mathbf {W}, \ beta) $ a aussi une distribution gaussienne, et $ p (\ mathbf {x}, \ mathbf {z} \ mid \ mathbf {W }, \ beta) $, qui est la marginalisation de la variable latente $ \ mathbf {z} $, $ p (\ mathbf {x} \ mid \ mathbf {W}, \ beta) $ a aussi une distribution gaussienne. En fait, c'est $ p (\ mathbf {x} \ mid \ mathbf {W}, \ beta) = \ mathcal {N} (\ mathbf {x} \ mid \ mathbb {E} [\ mathbf {x}], \ mathbb {V} [\ mathbf {x}]) $ peut être calculé comme suit.

\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}

Dans la transformation de formule ci-dessus, on suppose que $ \ mathbf {W} $ n'est pas une variable stochastique et peut être sorti de la valeur attendue, et que $ \ mathbf {z} $ et $ \ boldsymbol {\ epsilon} $ sont indépendants. J'utilise ça.

Ici, $ p (\ mathbf {x} \ mid \ mathbf {W}, \ beta) $ et $ \ mathbf { Il a la même distribution indépendante pour x} $. Donc

\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}

Il devient. Maximisez le logarithme de la fonction de vraisemblance dans l'équation (2) pour $ \ mathbf {W} $ et $ \ beta $ comme suit.

\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}

C'est la fonction objective de l'analyse probabiliste en composantes principales [^ 6]. Il y a deux façons d'estimer $ \ mathbf {W} $ et $ \ beta $ qui maximisent cette fonction objectif: en différenciant $ L $ et en la trouvant sous forme fermée, ou en utilisant l'algorithme EM. Puisqu'il est hors du sujet principal, il est omis ici. Veuillez consulter ce document pour plus de détails [^ 7]

Dual probabilistic PCA

Dans pPCA, nous avons considéré un modèle de probabilité de $ \ mathbf x_n $ avec $ \ mathbf z_n $ comme variable stochastique et $ \ mathbf {W} $ comme paramètre. Par contre, en PCA double probabiliste (DPPCA), le paramètre est $ \ mathbf {Z} = (\ mathbf z_1, \ mathbf z_2, \ cdots, \ mathbf z_N) \ in \ mathbb {R} ^ {L \ times N} $ ， $ \ Mathbf w_d $ est considéré comme une variable de probabilité Considérons un modèle probabiliste pour $ \ mathbf x_ {: d} $. Ici, $ \ mathbf x_ {: d} = (x_ {1d}, x_ {2d}, \ cdots, x_ {Nd}) ^ {\ rm T} $. En d'autres termes

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

Supposons que $ \ mathbf w_d $ et $ \ boldsymbol \ epsilon_ {: d} $ suivent chacun les distributions de probabilité suivantes.

\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) $ et $ p (\ mathbf x_ {: d} \ mid \ mathbf {Z}, \ mathbf {w} _ {: d}, \ beta) $ sont tous deux des distributions gaussiennes, le reste est donc pPCA De même, la fonction de vraisemblance peut être obtenue en marginalisant $ \ mathbf {w} _d $ comme suit.

\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}

La fonction objective de DPPCA peut être obtenue en prenant la valeur logarithmique de cette fonction de vraisemblance.

\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}

Ici, $ \ mathbf S $ et $ \ mathbf K $ sont définis comme suit.

\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}

En regardant l'équation (4), il peut sembler que ce que vous faites avec pPCA n'est pas si différent, car les rôles de $ \ mathbf W $ et $ \ mathbf Z $ sont juste échangés après tout. Cependant, considérons l'équation constante moins suivante (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}

A ce moment, l'estimation de $ \ mathbf {Z} $ qui maximise l'équation (4) à partir de l'équation suivante minimise la divergence KL entre la matrice gramme des données observées et la matrice gramme de la variable latente. Vous pouvez voir que cela équivaut à estimer {Z} $. En d'autres termes, le but de DPPCA est d'estimer la variable latente $ \ mathbf {Z} $ afin que la similitude entre les données observées et la similitude entre les variables latentes correspondantes correspondent autant que possible.

\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}

À ce stade, DPPCA peut réaliser une réduction de dimension non linéaire en définissant la similitude entre les données observées et la variable latente en utilisant la fonction noyau $ k (\ cdot, \ cdot) $ pour définir la similitude autre que le produit interne standard. .. En particulier, la méthode de réduction de dimension non linéaire à l'aide de la fonction noyau pour la similitude des données observées est appelée analyse en composantes principales du noyau [^ 8], et la méthode de réduction de dimension non linéaire utilisant la fonction noyau pour la similitude des variables latentes. Il s'appelle GPLVM.

Les deux méthodes permettent une réduction de dimension non linéaire, mais chacune a ses avantages et ses inconvénients. Puisque la fonction noyau est utilisée pour la similitude entre les données observées dans l'analyse en composantes principales du noyau, une fois que $ \ mathbf {S} $ est calculé en utilisant la fonction noyau, la solution analytique peut être obtenue de la même manière que l'analyse en composantes principales normales. .. Cependant, comme nous ne connaissons pas la cartographie de l'espace latent vers l'espace d'observation, il est nécessaire de résoudre le problème de pré-image [^ 9] afin d'estimer les points sur l'espace d'observation qui correspondent à des points quelconques de l'espace latent. D'autre part, GPLVM ne pose pas de problème de pré-image car il peut explicitement dessiner un mappage de l'espace latent vers l'espace d'observation. Au lieu de cela, vous devez mettre à jour $ \ mathbf {K} $ chaque fois que la variable latente change.

Gaussian Process Latent Variable Model

Jusqu'à présent, nous avons constaté que GPLVM est un modèle d'apprentissage qui estime la variable latente $ \ mathbf Z $ afin de maximiser la fonction objective suivante.

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}

Où $ \ mathbf K $ utilise la fonction du noyau $ k (\ cdot, \ cdot) $ vers $ \ mathbf K = k (\ mathbf Z, \ mathbf Z) + \ beta ^ {-1} \ mathbf {I Il est défini comme} _N $.

En fait, cette fonction objective peut également être dérivée du point de vue du médecin généraliste. Je vais omettre les détails sur GP cette fois, donc si vous voulez en savoir plus sur GP, veuillez vous référer à ce livre [^ 10]. Soit $ f: \ mathcal Z \ rightarrow \ mathcal X = \ mathbb {R} ^ D $ le mappage de l'espace latent $ \ mathcal Z $ vers l'espace d'observation $ \ mathcal X $, et $ \ mathbf x \ in \ mathcal X $ Est supposé être indépendant du point de vue dimensionnel. C'est,

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

On suppose que $ x_ {nd} $ est généré indépendamment pour chaque dimension à partir de $ \ mathbf z_n $ selon. Ici, $ \ epsilon_ {nd} $ est un bruit gaussien avec le paramètre de précision $ \ beta $. Soit la distribution antérieure de $ f_d $ $ f_d \ sim \ mathcal {GP} (0, k (\ mathbf {z}, \ mathbf {z} ')) $ [^ 11]. Lorsque les données observées sont $ \ mathbf X = (\ mathbf x_1, \ mathbf x_2, \ cdots, \ mathbf x_N) $ et la variable latente correspondante est $ \ mathbf Z $, la distribution précédente est $ p (\ mathbf f_d ). mid \ mathbf Z) = \ mathcal {N} (\ mathbf f_d \ mid \ mathbf 0, k (\ mathbf Z, \ mathbf Z)) $ [^ 12]. Ici, $ \ mathbf f_d = (f_d (\ mathbf z_1), f_d (\ mathbf z_2), \ cdots, f_d (\ mathbf z_N)) $. A ce moment, $ \ mathbf f_d $ est indépendant pour la dimension $ d $, donc la vraisemblance périphérique est

\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}

Il devient. De plus, $ p (\ mathbf {x} _ {: d} \ mid \ mathbf {f} _d, \ beta) $ et $ p (\ mathbf {f} _d \ mid \ mathbf {Z}) $ sont des distributions gaussiennes. Par conséquent, la vraisemblance périphérique a également une distribution gaussienne.

Donc $ p (\ mathbf x_ {: d} \ mid \ mathbf Z) = \ mathcal N (\ mathbf x_ {: d} \ mid \ mathbb E [\ mathbf x_ {: d}], \ mathbb V [\ mathbf x_ {: d}]) Si c'est $, ce sera comme suit.

\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}

Que ça

\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}

Et est d'accord avec l'équation (5). En d'autres termes, d'un point de vue différent, la variable latente $ \ mathbf Z $ estimée par DPPCA estime la variable latente qui maximise la vraisemblance marginale lorsqu'on considère un processus gaussien dans lequel la sortie est multidimensionnelle au lieu d'être une variable latente. Vous pouvez voir que c'est équivalent à cela. A partir de là, la distribution a posteriori du mapping $ f $ de l'espace latent vers l'espace d'observation peut s'écrire comme le processus gaussien suivant

\begin{align}
f &\sim \mathcal{GP}(\mu(\mathbf{z}),\sigma(\mathbf{z},\mathbf{z}')) \\
\mu(\mathbf{z}) &= k(\mathbf{z},\mathbf{Z}) (k(\mathbf{Z},\mathbf{Z})+\beta^{-1}\mathbf{I})^{-1}\mathbf{X} \\
\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}

De plus, concernant les hyperparamètres, on peut voir qu'estimer les hyperparamètres qui maximisent l'équation (5) équivaut à estimer les hyperparamètres qui maximisent la vraisemblance périphérique, donc la validité de l'estimation des hyperparamètres dans l'apprentissage non supervisé La sexualité peut également être garantie dans l'optique de maximiser la probabilité périphérique.

Puisque l'équation (5) ne peut pas être résolue analytiquement, $ \ mathbf Z $ et les hyperparamètres sont estimés en utilisant la méthode du gradient. Lors de la différenciation de $ L_ {\ rm DPPCA} $ pour $ z_ {nl} $, $ z_ {nl} $ est un argument de la fonction noyau, il peut donc être calculé en utilisant la loi de chaîne différentielle suivante. Je peux le faire.

\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}

ici

\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}

Il devient. Pour $ \ frac {\ partial \ mathbf {K}} {\ partial z_ {nl}} $, cela dépend de la fonction du noyau, donc différenciez en conséquence. De plus, dans le calcul de l'algorithme actuel, la formule avec le logarithmique $ \ log {p (\ mathbf {Z})} $ de la distribution précédente est maximisée afin d'empêcher la variable latente de devenir une valeur extrême. La distribution normale standard est essentiellement utilisée comme distribution antérieure.

Vérification du fonctionnement de GPLVM

Enfin, nous allons réellement implémenter GPLVM sur python et vérifier si l'apprentissage fonctionne avec des données simples. Cette fois, les données obtenues en échantillonnant au hasard 200 points de la variable latente $ \ mathbf z $ dans la plage de $ [-1,1] ^ 2 $ ont été mappées à l'espace tridimensionnel par la fonction suivante en tant que données d'observation. Où $ \ mathbf {\ epsilon} $ est le bruit gaussien. Le noyau gaussien est utilisé comme fonction de noyau, et les hyperparamètres et le bruit observé sont estimés en maximisant la vraisemblance périphérique. La valeur initiale de la variable latente est également déterminée aléatoirement.

\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}

Le programme actuel est le suivant.

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):

            #Mise à jour des variables latentes
            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 = -(((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)
            self.X = self.X + epsilonX * dFdX
            self.history['X'][i] = self.X

            #Mise à jour des hyper paramètres
            Dist = ((self.X[:, None, :] - self.X[None, :, :]) ** 2).sum(axis=2)
            dKdSigma = 0.5*Dist/self.hyperparam[0]*K
            dFdSigma = np.trace(G @ dKdSigma)
            sigma = sigma + epsilonSigma * dFdSigma
            self.hyperparam[0] = np.exp(sigma)

            dKdAlpha = self.hyperparam[1]*np.eye(self.dataNum)
            dFdAlpha = np.trace(G @ dKdAlpha)
            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

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, isSave=False):
    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', aspect='equal')
            ax_input.cla()
            #Affichage de l'espace d'observation
            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())
            # plt.savefig("fig1.pdf")

            #Affichage de l'espace latent
            ax_latent = fig.add_subplot(1, 2, 2, aspect='equal')
            ax_latent.cla()
            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("fig/fig{0}.png ".format(i))

        plt.pause(0.001)

    if isSave:
        plt.savefig("result.png ", dpi=100)
    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

    #Génération de données d'entrée
    # np.random.seed(seedData)
    [X,Y] = createKuraData(N,D,sigma=0.01)
    # Y = np.loadtxt('kura.txt', delimiter=' ')

    #Paramètres du noyau
    [U,D,Vt] = np.linalg.svd(Y)
    model = GPLVM(Y,L, np.array([sigma**2,alpha/beta]))
    #Optimisation GPLVM
    epoch = 200
    model.fit(epoch=epoch,epsilonX=0.05,epsilonSigma=0.0005,epsilonAlpha=0.00001)

    #Obtenez la variable latente estimée
    X = model.history['X']
    f = model.history['F']

    #Affichage des résultats d'apprentissage
    plot_prediction(Y,f,X,epoch,True)

La figure de gauche montre le processus d'apprentissage de la diversité estimée dans l'espace d'observation. Le chiffre de droite est la variable latente à ce moment. Si la valeur initiale du paramètre n'est pas mauvaise dans une certaine mesure, il fonctionne de manière stable même si la valeur initiale est déterminée de manière aléatoire. En particulier, il a tendance à être plus facile à stabiliser si la largeur de la fonction de noyau est fixée à l'avance.

en conclusion

J'ai résumé ma compréhension de GPLVM. Nous vous serions reconnaissants si vous pouviez nous contacter si vous avez des questions ou des suggestions.