Letztes Mal wurde geändert. Der Inhalt ist wie unten angegeben.
Reduzierung der Speichernutzung
Einführung einer Fehlerbehandlung
Die Beschleunigung wird später kommen.
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.
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.
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.
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.
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.
»Es ist länger als ich erwartet hatte.
――Die Funktion einsum ist ebenfalls nützlich.
Recommended Posts