[PYTHON] Ich habe versucht, die Grundform von GPLVM zusammenzufassen

Dieser Artikel stammt aus Furukawa Lab Advent_calendar Day 23.

Einführung

GPLVM ist eine unbeaufsichtigte Lernmethode unter Verwendung des Gaußschen Prozesses (GP), mit der niedrigdimensionale Varianten geschätzt werden können, bei denen Daten im hochdimensionalen Raum verteilt werden. Der attraktivste Punkt von GPLVM ist, dass das Modell sehr einfach ist, daher sehr erweiterbar und theoretisch einfach zu handhaben ist. Tatsächlich wurde es auf komplexere Situationen wie die Schätzung von Hyperparametern im Rahmen der Bayes'schen Schätzung, der Zeitreihendatenanalyse [^ 1], der Multiview-Analyse [^ 2] und der Tensorzerlegung [^ 3] ausgedehnt. .. Diesmal wird der grundlegendste GPLVM-Algorithmus ausgehend von der probabilistischen PCA (pPCA) abgeleitet. Diese Ableitung folgt im Allgemeinen dem Ablauf der Arbeit von Professor Lawrence. Bitte wenden Sie sich für Einzelheiten an diese Person [^ 4].

probabilistic PCA

Wie der Name schon sagt, ist pPCA ein Modell, das die Hauptkomponentenanalyse im Rahmen der Wahrscheinlichkeitstheorie neu verteilt. In pPCA ist $ \ mathbf {X} = (\ mathbf {x} _1, \ mathbf {x} _2, \ cdots, \ mathbf {x} _N) \ in \ mathbb {R} ^ {D \ times N} $ Angenommen, $ \ mathbf {x} $ wird beobachtet, indem die latente Variable $ \ mathbf {z} \ in \ mathbb {R} ^ L $ wie folgt zugeordnet wird [^ 5] ].

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

Hier ist $ \ epsilon $ das beobachtete Rauschen, und das isotrope Gaußsche Rauschen mit dem Genauigkeitsparameter $ \ beta ^ {-1} $ wird angenommen. Es ist auch $ \ mathbf {W} \ in \ mathbb {R} ^ {D \ times L} $. Wir nehmen auch an, dass $ \ mathbf {z} $ einer Standardnormalverteilung der $ L $ -Dimension folgt. Das heißt, $ \ mathbf {x} $ und $ \ mathbf {z} $ folgen der folgenden Wahrscheinlichkeitsverteilung.

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

Zu diesem Zeitpunkt sind sowohl $ p (\ mathbf {x} \ mid \ mathbf {W}, \ mathbf {z}, \ beta) $ als auch $ p (\ mathbf {z}) $ Gaußsche Verteilungen, also gleichzeitige Verteilung $ p (\ mathbf {x}, \ mathbf {z} \ mid \ mathbf {W}, \ beta) $ hat ebenfalls eine Gaußsche Verteilung und $ p (\ mathbf {x}, \ mathbf {z} \ mid \ mathbf {W. }, \ beta) $, die Marginalisierung der latenten Variablen $ \ mathbf {z} $, $ p (\ mathbf {x} \ mid \ mathbf {W}, \ beta) $, hat ebenfalls eine Gaußsche Verteilung. Tatsächlich ist dies $ p (\ mathbf {x} \ mid \ mathbf {W}, \ beta) = \ mathcal {N} (\ mathbf {x} \ mid \ mathbb {E} [\ mathbf {x}], \ mathbb {V} [\ mathbf {x}]) $ kann wie folgt berechnet werden.

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

In der obigen Formeltransformation wird angenommen, dass $ \ mathbf {W} $ keine stochastische Variable ist und aus dem erwarteten Wert herausgesetzt werden kann und dass $ \ mathbf {z} $ und $ \ boldsymbol {\ epsilon} $ unabhängig sind. Ich benutze das.

Hier $ p (\ mathbf {x} \ mid \ mathbf {W}, \ beta) $ und $ \ mathbf { Es hat die gleiche unabhängige Verteilung für x} $. Deshalb

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

Es wird. Maximieren Sie den Logarithmus der Wahrscheinlichkeitsfunktion in Gleichung (2) für $ \ mathbf {W} $ und $ \ beta $ wie folgt.

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

Dies ist die objektive Funktion der probabilistischen Hauptkomponentenanalyse [^ 6]. Es gibt zwei Möglichkeiten, $ \ mathbf {W} $ und $ \ beta $ zu schätzen, die diese Zielfunktion maximieren: durch Differenzieren von $ L $ und Finden in geschlossener Form oder Verwenden des EM-Algorithmus. Da es nicht zum Hauptthema gehört, wird es hier weggelassen. Weitere Informationen finden Sie in diesem Dokument. [^ 7]

Dual probabilistic PCA

In pPCA haben wir ein Wahrscheinlichkeitsmodell von $ \ mathbf x_n $ mit $ \ mathbf z_n $ als stochastischer Variable und $ \ mathbf {W} $ als Parameter betrachtet. Andererseits ist in Dual Probabilistic PCA (DPPCA) der Parameter $ \ mathbf {Z} = (\ mathbf z_1, \ mathbf z_2, \ cdots, \ mathbf z_N) \ in \ mathbb {R} ^ {L \ times N} $ ， $ \ Mathbf w_d $ wird als Wahrscheinlichkeitsvariable angesehen Betrachten Sie ein Wahrscheinlichkeitsmodell für $ \ mathbf x_ {: d} $. Hier ist $ \ mathbf x_ {: d} = (x_ {1d}, x_ {2d}, \ cdots, x_ {Nd}) ^ {\ rm T} $. Mit anderen Worten

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

Nehmen wir an, dass $ \ mathbf w_d $ und $ \ boldsymbol \ epsilon_ {: d} $ jeweils den folgenden Wahrscheinlichkeitsverteilungen folgen.

\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) $ und $ p (\ mathbf x_ {: d} \ mid \ mathbf {Z}, \ mathbf {w} _ {: d}, \ beta) $ sind beide Gaußsche Verteilungen In ähnlicher Weise kann die Wahrscheinlichkeitsfunktion erhalten werden, indem $ \ mathbf {w} _d $ wie folgt marginalisiert wird.

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

Die Zielfunktion von DPPCA kann erhalten werden, indem der logarithmische Wert dieser Wahrscheinlichkeitsfunktion genommen wird.

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

Hier sind $ \ mathbf S $ und $ \ mathbf K $ wie folgt definiert.

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

Wenn man sich Gleichung (4) ansieht, scheint es, dass das, was Sie mit pPCA machen, nicht so unterschiedlich ist, da die Rollen von $ \ mathbf W $ und $ \ mathbf Z $ doch nur vertauscht werden. Betrachten Sie jedoch die folgende Konstante minus Gleichung (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}

Zu diesem Zeitpunkt minimiert das Schätzen von $ \ mathbf {Z} $, das Gleichung (4) aus der folgenden Gleichung maximiert, die KL-Divergenz zwischen der Grammmatrix der beobachteten Daten und der Grammmatrix der latenten Variablen $ \ mathbf. Sie können sehen, dass dies der Schätzung von {Z} $ entspricht. Mit anderen Worten, der Zweck von DPPCA besteht darin, die latente Variable $ \ mathbf {Z} $ so zu schätzen, dass die Ähnlichkeit zwischen den beobachteten Daten und die Ähnlichkeit zwischen den entsprechenden latenten Variablen so weit wie möglich übereinstimmen.

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

Zu diesem Zeitpunkt kann DPPCA eine nichtlineare Dimensionsreduktion realisieren, indem die Ähnlichkeit der beobachteten Daten und die Ähnlichkeit latenter Variablen unter Verwendung der Kernelfunktion $ k (\ cdot, \ cdot) $ definiert werden, um die Ähnlichkeit außer dem inneren Standardprodukt zu definieren. .. Insbesondere wird das Verfahren zum Durchführen einer nichtlinearen Dimensionsreduktion unter Verwendung der Kernelfunktion für die Ähnlichkeit der beobachteten Daten als Kernel-Hauptkomponentenanalyse [^ 8] bezeichnet, und das Verfahren zum Durchführen einer nichtlinearen Dimensionsreduktion unter Verwendung der Kernelfunktion für die Ähnlichkeit latenter Variablen wird aufgerufen. Es heißt GPLVM.

Beide Methoden ermöglichen eine nichtlineare Dimensionsreduzierung, haben jedoch jeweils ihre Vor- und Nachteile. Da die Kernelfunktion für die Ähnlichkeit zwischen den beobachteten Daten in der Kernel-Hauptkomponentenanalyse verwendet wird, kann die analytische Lösung auf die gleiche Weise wie die normale Hauptkomponentenanalyse erhalten werden, sobald $ \ mathbf {S} $ unter Verwendung der Kernelfunktion berechnet wurde. .. Da wir jedoch die Abbildung vom latenten Raum auf den Beobachtungsraum nicht kennen, ist es notwendig, das Vorbildproblem [^ 9] zu lösen, um die Punkte auf dem Beobachtungsraum zu schätzen, die irgendwelchen Punkten auf dem latenten Raum entsprechen. Andererseits verursacht GPLVM das Vorbildproblem nicht, da es explizit eine Zuordnung vom latenten Raum zum Beobachtungsraum zeichnen kann. Stattdessen müssen Sie $ \ mathbf {K} $ jedes Mal aktualisieren, wenn sich die latente Variable ändert.

Gaussian Process Latent Variable Model

Bisher haben wir festgestellt, dass GPLVM ein Lernmodell ist, das die latente Variable $ \ mathbf Z $ schätzt, um die folgende Zielfunktion zu maximieren.

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}

Wobei $ \ mathbf K $ die Kernelfunktion $ k (\ cdot, \ cdot) $ zu $ \ mathbf K = k (\ mathbf Z, \ mathbf Z) + \ beta ^ {-1} \ mathbf {I verwendet Es ist definiert als} _N $.

Tatsächlich kann diese Zielfunktion auch aus der Sicht von GP abgeleitet werden. Ich werde diesmal die Details zu GP weglassen. Wenn Sie also etwas über GP wissen möchten, lesen Sie bitte dieses Buch [^ 10]. Die Zuordnung vom latenten Raum $ \ mathcal Z $ zum Beobachtungsraum $ \ mathcal X $ sei $ f: \ mathcal Z \ rightarrow \ mathcal X = \ mathbb {R} ^ D $ und $ \ mathbf x \ in \ mathcal X. $ Wird als dimensionsunabhängig angenommen. Das ist,

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

Es wird angenommen, dass $ x_ {nd} $ unabhängig für jede Dimension aus $ \ mathbf z_n $ gemäß generiert wird. Hier ist $ \ epsilon_ {nd} $ Gaußsches Rauschen mit dem Präzisionsparameter $ \ beta $. Die vorherige Verteilung von $ f_d $ sei $ f_d \ sim \ mathcal {GP} (0, k (\ mathbf {z}, \ mathbf {z} ')) $ [^ 11]. Wenn die beobachteten Daten $ \ mathbf X = (\ mathbf x_1, \ mathbf x_2, \ cdots, \ mathbf x_N) $ sind und die entsprechende latente Variable $ \ mathbf Z $ ist, ist die vorherige Verteilung $ p (\ mathbf f_d ). mid \ mathbf Z) = \ mathcal {N} (\ mathbf f_d \ mid \ mathbf 0, k (\ mathbf Z, \ mathbf Z)) $ [^ 12]. Hier $ $ mathbf f_d = (f_d (\ mathbf z_1), f_d (\ mathbf z_2), \ cdots, f_d (\ mathbf z_N)) $. Zu diesem Zeitpunkt ist $ \ mathbf f_d $ unabhängig von der Dimension $ d $, also ist die periphere Wahrscheinlichkeit

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

Es wird. Außerdem sind $ p (\ mathbf {x} _ {: d} \ mid \ mathbf {f} _d, \ beta) $ und $ p (\ mathbf {f} _d \ mid \ mathbf {Z}) $ Gaußsche Verteilungen. Daher hat die periphere Wahrscheinlichkeit auch eine Gaußsche Verteilung.

Daher ist $ p (\ mathbf x_ {: d} \ mid \ mathbf Z) = \ mathcal N (\ mathbf x_ {: d} \ mid \ mathbb E [\ mathbf x_ {: d}], \ mathbb V [\ mathbf x_ {: d}]) Wenn es $ ist, ist es wie folgt.

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

Als das

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

Und stimmt mit Gleichung (5) überein. Mit anderen Worten, aus einem anderen Blickwinkel schätzt die von DPPCA geschätzte latente Variable $ \ mathbf Z $ die latente Variable, die die Grenzwahrscheinlichkeit maximiert, wenn ein Gaußscher Prozess betrachtet wird, bei dem die Ausgabe mehrdimensional ist, anstatt eine latente Variable zu sein. Sie können sehen, dass es dem entspricht. Daraus kann die posteriore Verteilung der Abbildung $ f $ vom latenten Raum zum Beobachtungsraum als der folgende Gaußsche Prozess geschrieben werden.

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

Darüber hinaus ist in Bezug auf Hyperparameter ersichtlich, dass das Schätzen der Hyperparameter, die Gleichung (5) maximieren, dem Schätzen der Hyperparameter entspricht, die die periphere Wahrscheinlichkeit maximieren, so dass die Hyperparameterschätzung beim unbeaufsichtigten Lernen angemessen ist. Sexualität kann auch unter dem Gesichtspunkt der Maximierung der peripheren Wahrscheinlichkeit garantiert werden.

Da Gleichung (5) nicht analytisch gelöst werden kann, werden $ \ mathbf Z $ und Hyperparameter unter Verwendung der Gradientenmethode geschätzt. Bei der Differenzierung von $ L_ {\ rm DPPCA} $ nach $ z_ {nl} $ ist $ z_ {nl} $ ein Argument der Kernelfunktion, sodass es unter Verwendung des folgenden Differentialkettengesetzes berechnet werden kann. Ich kann es schaffen

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

Hier

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

Es wird. Für $ \ frac {\ partielle \ mathbf {K}} {\ partielle z_ {nl}} $ hängt es von der Kernelfunktion ab, also differenzieren Sie entsprechend. Außerdem wird bei der tatsächlichen Algorithmusberechnung die Formel mit dem logarithmischen $ \ log {p (\ mathbf {Z})} $ der vorherigen Verteilung maximiert, um zu verhindern, dass die latente Variable zu einem Extremwert wird. Die Standardnormalverteilung wird grundsätzlich als vorherige Verteilung verwendet.

Betriebsüberprüfung von GPLVM

Schließlich werden wir GPLVM tatsächlich auf Python implementieren und überprüfen, ob das Lernen mit einfachen Daten funktioniert. Dieses Mal wurden die Daten, die durch zufälliges Abtasten von 200 Punkten der latenten Variablen $ \ mathbf z $ im Bereich von $ [-1,1] ^ 2 $ erhalten wurden, durch die folgende Funktion als Beobachtungsdaten auf den dreidimensionalen Raum abgebildet. Wobei $ \ mathbf {\ epsilon} $ Gaußsches Rauschen ist. Der Gaußsche Kernel wird als Kernelfunktion verwendet, und Hyperparameter und beobachtetes Rauschen werden durch Maximieren der peripheren Wahrscheinlichkeit geschätzt. Der Anfangswert der latenten Variablen wird ebenfalls zufällig bestimmt.

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

Das eigentliche Programm ist wie folgt.

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

            #Aktualisierung der latenten Variablen
            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

            #Hyper-Parameter-Update
            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()
            #Anzeige des Beobachtungsraums
            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")

            #Anzeige des latenten Raums
            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

    #Generierung von Eingabedaten
    # np.random.seed(seedData)
    [X,Y] = createKuraData(N,D,sigma=0.01)
    # Y = np.loadtxt('kura.txt', delimiter=' ')

    #Kernel-Einstellungen
    [U,D,Vt] = np.linalg.svd(Y)
    model = GPLVM(Y,L, np.array([sigma**2,alpha/beta]))
    #GPLVM-Optimierung
    epoch = 200
    model.fit(epoch=epoch,epsilonX=0.05,epsilonSigma=0.0005,epsilonAlpha=0.00001)

    #Holen Sie sich die geschätzte latente Variable
    X = model.history['X']
    f = model.history['F']

    #Anzeige der Lernergebnisse
    plot_prediction(Y,f,X,epoch,True)

Die Abbildung links zeigt den Lernprozess der geschätzten Diversität im Beobachtungsraum. Die Abbildung rechts ist die latente Variable zu diesem Zeitpunkt. Wenn der Anfangswert des Parameters bis zu einem gewissen Grad nicht schlecht ist, arbeitet er stabil, selbst wenn der Anfangswert zufällig bestimmt wird. Insbesondere ist es tendenziell einfacher zu stabilisieren, wenn die Breite der Kernelfunktion im Voraus groß eingestellt wird.

abschließend

Ich habe mein Verständnis von GPLVM zusammengefasst. Wir würden uns freuen, wenn Sie uns kontaktieren könnten, wenn Sie Fragen oder Anregungen haben.