Ich habe versucht, PLSA in Python 2 zu implementieren

Überblick

Letztes Mal wurde geändert. Der Inhalt ist wie unten angegeben.

Die Beschleunigung wird später kommen.

Reduzierung der Speichernutzung

Ursache und Politik

Erstens ist die vorherige Implementierung sehr speicherintensiv. Die Ursache ist, dass P (z | x, y) im E-Schritt des EM-Algorithmus berechnet wird. Da wir einfach ein dreidimensionales Array erstellt haben, wird die Speichernutzung erheblich reduziert, wenn wir es nicht erstellen.

Schauen wir uns insbesondere die Aktualisierungsformel für P (x | z) im Schritt M an. Da der Nenner nur normalisiert ist, betrachten wir nur das Molekül.

P\left( x | z \right)Molekül= \sum_{y} N_{x, y} P \left( z | x, y \right)

Ersetzen Sie die Aktualisierungsformel für P (z | x, y) in Schritt E durch diese.

\begin{equation}
P\left( x | z \right)Molekül= \sum_{y} N_{x, y} \frac{P\left( z \right)P\left( x | z \right)P\left( y | z \right)}{\sum_{z} P\left( z \right)P\left( x | z \right)P\left( y | z \right)} \tag{1}
\end{equation}

Ich werde diesen Ausdruck implementieren, aber ich möchte die for-Anweisung nicht herumdrehen, also verwende ich numpys einsum.

numpy.einsum()

Die Einsum-Funktion ist eine Reduktion von Einstein. Es ist schwer zu verstehen, deshalb hier ein Beispiel.

P(x,y) = \sum_{z}P(z)P(x|z)P(y|z)

Wenn Sie implementieren

Pxy = numpy.einsum('k,ki,kj->ij', Pz, Px_z, Py_z)

Es wird sein.

Gleichung (1) wird unter Verwendung dieser einsum-Funktion implementiert, aber es ist schwierig, sie so zu implementieren, wie sie ist, so dass die Gleichung wie folgt transformiert wird.

P\left( x | z \right)Molekül= \sum_{y} \frac{N_{x, y}}{\sum_{z} P\left( z \right)P\left( x | z \right)P\left( y | z \right)} P\left( z \right)P\left( x | z \right)P\left( y | z \right)

Wenn Sie dies implementieren, sieht es so aus.

tmp = N / numpu.einsum('k,ki,kj->ij', Pz, Px_z, Py_z)
Px_z = numpy.einsum('ij,k,ki,kj->ki', tmp, Pz, Px_z, Py_z)

Wir werden später sehen, wie viel Speicherplatz dadurch reduziert wurde.

Fehlerbehandlung

Ich habe es mit der Funktion einsum implementiert,

tmp = N / numpu.einsum('k,ki,kj->ij', Pz, Px_z, Py_z)

Berücksichtigen Sie die Fehlerbehandlung, wenn Sie durch 0 in geteilt werden. Die Tatsache, dass dieser Nenner 0 ist, bedeutet, dass für einige x und y,

\sum_{z}P(z)P(x|z)P(y|z) = 0

Es kommt also kein negativer Wert heraus, also für einige x und y und alle z,

P(z)P(x|z)P(y|z) = 0

Ist festgelegt. Mit anderen Worten ist der E-Schritt des EM-Algorithmus wie folgt.

\begin{eqnarray}
P(z|x,y) & = & \frac{P\left( z \right)P\left( x | z \right)P\left( y | z \right)}{\sum_{z} P\left( z \right)P\left( x | z \right)P\left( y | z \right)}
& = & 0
\end{eqnarray}

Daher wird das durch 0 geteilte Element auf 0 gesetzt.

Wo in numpy

1 / 0 = inf
0 / 0 = nan

Verwenden Sie also numpy.isinf bzw. numpy.isnan

tmp = N / numpu.einsum('k,ki,kj->ij', Pz, Px_z, Py_z)
tmp[numpy.isinf(tmp)] = 0
tmp[numpy.isnan(tmp)] = 0

Px_z = numpy.einsum('ij,k,ki,kj->ki', tmp, Pz, Px_z, Py_z)

Es wird sein.

Implementierung

Zusammenfassend ist die Gesamtimplementierung wie folgt.

plsa.py


import numpy as np


class PLSA(object):
    def __init__(self, N, Z):
        self.N = N
        self.X = N.shape[0]
        self.Y = N.shape[1]
        self.Z = Z

        # P(z)
        self.Pz = np.random.rand(self.Z)
        # P(x|z)
        self.Px_z = np.random.rand(self.Z, self.X)
        # P(y|z)
        self.Py_z = np.random.rand(self.Z, self.Y)

        #Normalisierung
        self.Pz /= np.sum(self.Pz)
        self.Px_z /= np.sum(self.Px_z, axis=1)[:, None]
        self.Py_z /= np.sum(self.Py_z, axis=1)[:, None]

    def train(self, k=200, t=1.0e-7):
        '''
Wiederholen Sie die Schritte E und M, bis die Protokollwahrscheinlichkeit konvergiert
        '''
        prev_llh = 100000
        for i in xrange(k):
            self.em_algorithm()
            llh = self.llh()

            if abs((llh - prev_llh) / prev_llh) < t:
                break

            prev_llh = llh

    def em_algorithm(self):
        '''
EM-Algorithmus
        P(z), P(x|z), P(y|z)Aktualisieren
        '''
        tmp = self.N / np.einsum('k,ki,kj->ij', self.Pz, self.Px_z, self.Py_z)
        tmp[np.isnan(tmp)] = 0
        tmp[np.isinf(tmp)] = 0

        Pz = np.einsum('ij,k,ki,kj->k', tmp, self.Pz, self.Px_z, self.Py_z)
        Px_z = np.einsum('ij,k,ki,kj->ki', tmp, self.Pz, self.Px_z, self.Py_z)
        Py_z = np.einsum('ij,k,ki,kj->kj', tmp, self.Pz, self.Px_z, self.Py_z)

        self.Pz = Pz / np.sum(Pz)
        self.Px_z = Px_z / np.sum(Px_z, axis=1)[:, None]
        self.Py_z = Py_z / np.sum(Py_z, axis=1)[:, None]

    def llh(self):
        '''
Log-Wahrscheinlichkeit
        '''
        Pxy = np.einsum('k,ki,kj->ij', self.Pz, self.Px_z, self.Py_z)
        Pxy /= np.sum(Pxy)
        lPxy = np.log(Pxy)
        lPxy[np.isinf(lPxy)] = -1000

        return np.sum(self.N * lPxy)

Bei der Berechnung der Protokollwahrscheinlichkeit kann ein Unterlauf auftreten, der zu log (0) = -inf führt. Da der Mindestwert der Gleitkommazahl mit doppelter Genauigkeit etwa 4,94e-324 beträgt, ist er kleiner als log (4,94e-324) = -744,4, sodass ungefähr -1000 eingegeben wird.

Messung

Verwenden Sie memory_profiler, um zu messen, wie viel Speicher belegt wurde. Ich habe es mit dem folgenden Skript gemessen.

memory_profile.py


import numpy as np
from memory_profiler import profile

X = 1000
Y = 1000
Z = 10


@profile
def main():
    from plsa import PLSA
    plsa = PLSA(np.random.rand(X, Y), Z)
    plsa.em_algorithm()
    llh = plsa.llh()


if __name__ == '__main__':
    main()

Im Fall von X = 1000, Y = 1000, Z = 10 in der vorherigen Implementierung,

$ python profile_memory_element_wise_product.py 
Filename: profile_memory_element_wise_product.py

Line #    Mem usage    Increment   Line Contents
================================================
    10     15.9 MiB      0.0 MiB   @profile
    11                             def main():
    12     15.9 MiB      0.0 MiB       from plsa_element_wise_product import PLSA
    13     23.9 MiB      8.0 MiB       plsa = PLSA(np.random.rand(X, Y), Z)
    14    108.0 MiB     84.1 MiB       plsa.e_step()
    15    184.5 MiB     76.5 MiB       plsa.m_step()
    16    199.8 MiB     15.3 MiB       llh = plsa.llh()

In dieser Implementierung

$ python profile_memory_einsum.py 
Filename: profile_memory_einsum.py

Line #    Mem usage    Increment   Line Contents
================================================
    10     15.8 MiB      0.0 MiB   @profile
    11                             def main():
    12     15.8 MiB      0.0 MiB       from plsa_einsum import PLSA
    13     23.7 MiB      7.9 MiB       plsa = PLSA(np.random.rand(X, Y), Z)
    14     40.7 MiB     16.9 MiB       plsa.em_algorithm()
    15     48.4 MiB      7.8 MiB       llh = plsa.llh()

Im Fall von X = 5000, Y = 5000, Z = 10 in der vorherigen Implementierung,

$ python profile_memory_element_wise_product.py 
Filename: profile_memory_element_wise_product.py

Line #    Mem usage    Increment   Line Contents
================================================
    10     15.9 MiB      0.0 MiB   @profile
    11                             def main():
    12     15.9 MiB      0.0 MiB       from plsa_element_wise_product import PLSA
    13    207.6 MiB    191.7 MiB       plsa = PLSA(np.random.rand(X, Y), Z)
    14   2115.4 MiB   1907.8 MiB       plsa.e_step()
    15   2115.5 MiB      0.1 MiB       plsa.m_step()
    16   2115.5 MiB      0.0 MiB       llh = plsa.llh()

In dieser Implementierung

$ python profile_memory_einsum.py 
Filename: profile_memory_einsum.py

Line #    Mem usage    Increment   Line Contents
================================================
    10     15.7 MiB      0.0 MiB   @profile
    11                             def main():
    12     15.7 MiB      0.0 MiB       from plsa_einsum import PLSA
    13    207.5 MiB    191.7 MiB       plsa = PLSA(np.random.rand(X, Y), Z)
    14    233.0 MiB     25.6 MiB       plsa.em_algorithm()
    15    233.1 MiB      0.0 MiB       llh = plsa.llh()

Zusammenfassend ist die gesamte Speichernutzung

Implementierung X=1000,Y=1000,Z=10 X=5000,Y=5000,Z=10
Letzte Implementierung 199.8 MiB 2115.5 MiB
Diese Implementierung 48.4 MiB 233.1 MiB

Wenn wir dies jedoch auf den EM-Algorithmus und den Berechnungsteil der logarithmischen Wahrscheinlichkeit beschränken,

Implementierung X=1000,Y=1000,Z=10 X=5000,Y=5000,Z=10
Letzte Implementierung 175.9 MiB 1907.9 MiB
Diese Implementierung 24.7 MiB 25.6 MiB

Es ist ersichtlich, dass die in dieser Implementierung verwendete Speichermenge kaum zugenommen hat. Selbst wenn die Anzahl der Daten zunimmt, habe ich damit keine Angst vor dem Speicher.

Berechnungsgeschwindigkeit

Die einsum-Funktion ist praktisch, aber es dauert lange, alle drei oder vier gleichzeitig wie diese zu berechnen. Bei meinem MacBook dauerte es ungefähr 8,7 Sekunden, um die Protokollwahrscheinlichkeit und den EM-Algorithmus einmal bei X = 5000, Y = 5000, Z = 10 zu berechnen. Dies ist nicht realistisch und muss verbessert werden, aber es wird zu einem späteren Zeitpunkt mehr geben.

Impressionen

»Es ist länger als ich erwartet hatte.

――Die Funktion einsum ist ebenfalls nützlich.

Recommended Posts

Ich habe versucht, PLSA in Python zu implementieren
Ich habe versucht, PLSA in Python 2 zu implementieren
Ich habe versucht, Permutation in Python zu implementieren
Ich habe versucht, ADALINE in Python zu implementieren
Ich habe versucht, PPO in Python zu implementieren
Ich habe versucht, TOPIC MODEL in Python zu implementieren
Ich habe versucht, eine selektive Sortierung in Python zu implementieren
Ich habe versucht, einen Pseudo-Pachislot in Python zu implementieren
Ich habe versucht, Drakues Poker in Python zu implementieren
Ich habe versucht, GA (genetischer Algorithmus) in Python zu implementieren
Ich habe versucht, einen eindimensionalen Zellautomaten in Python zu implementieren
Ich habe versucht, die Mail-Sendefunktion in Python zu implementieren
Ich habe versucht, das Blackjack of Trump-Spiel mit Python zu implementieren
Ich habe versucht, PCANet zu implementieren
Ich habe versucht, ein missverstandenes Gefangenendilemma in Python zu implementieren
Ich habe versucht, StarGAN (1) zu implementieren.
Ich habe versucht, die Bayes'sche lineare Regression durch Gibbs-Sampling in Python zu implementieren
Ich habe versucht, Trumps Kartenspiel in Python zu implementieren
Ich habe versucht, die in Python installierten Pakete grafisch darzustellen
Ich möchte Timeout einfach in Python implementieren
Ich habe versucht, Mine Sweeper auf dem Terminal mit Python zu implementieren
Ich habe versucht, künstliches Perzeptron mit Python zu implementieren
Ich habe versucht zusammenzufassen, wie man Pandas von Python benutzt
Ich habe versucht, Deep VQE zu implementieren
Ich habe versucht, Python zu berühren (Installation)
Ich habe versucht, eine kontroverse Validierung zu implementieren
Ich habe versucht, Realness GAN zu implementieren
Ich habe Line Benachrichtigung in Python versucht
Ich habe versucht, die Zusammenführungssortierung in Python mit möglichst wenigen Zeilen zu implementieren
Ich habe versucht, ein scheinbar Windows-Snipper-Tool mit Python zu implementieren
Ich habe versucht, API list.csv mit Python aus swagger.yaml zu erstellen
Ich habe versucht "Wie man eine Methode in Python dekoriert"
Ich habe eine Stoppuhr mit tkinter mit Python gemacht
Ich habe versucht, die Behandlung von Python-Ausnahmen zusammenzufassen
Ich habe versucht, Autoencoder mit TensorFlow zu implementieren
Python3-Standardeingabe habe ich versucht zusammenzufassen
Ich habe versucht, die Bayes'sche Optimierung von Python zu verwenden
Ich wollte ABC159 mit Python lösen
Ich habe versucht, CVAE mit PyTorch zu implementieren
[Python] Ich habe versucht, TF-IDF stetig zu berechnen
Ich habe versucht, Python zu berühren (grundlegende Syntax)
[Python] Ich habe versucht, eine stabile Sortierung zu implementieren
Ich habe Python> autopep8 ausprobiert
Implementieren Sie XENO mit Python
Ich habe versucht zu debuggen.
Implementieren Sie sum in Python
Ich habe Python> Decorator ausprobiert
Implementieren Sie Traceroute in Python 3
Ich habe versucht, das Lesen von Dataset mit PyTorch zu implementieren
Ich möchte Dunnetts Test in Python machen
Versuchen Sie, Oni Mai Tsuji Miserable mit Python zu implementieren
Python: Ich konnte in Lambda rekursieren
Ich möchte mit Python ein Fenster erstellen
Ich habe versucht, mit Python ein Tippspiel zu spielen
So implementieren Sie Shared Memory in Python (mmap.mmap)
Ich habe die Methode der kleinsten Quadrate in Python ausprobiert