Ajuste de un modelo de mezcla gaussiana usando descenso de gradiente estocástico

Estoy trabajando en un modelo de aprendizaje de categoría en línea que utiliza el descenso de gradiente estocástico para adaptarse a un modelo de mezcla gaussiana. El modelo se basa en el modelo de aprendizaje en línea utilizado en Toscano y McMurray (2010).

Si bien el descenso de gradiente parece estar funcionando bastante bien para estimar las medias y frecuencias / probabilidades de mezcla de las categorías, estoy teniendo problemas para estimar las covarianzas de los componentes de la mezcla. Las derivadas parciales que he estado usando para la actualización de descenso de gradiente provienen de Petersen y Pedersen (2008) (p. 44)

Empezando con

$p(x) = \sum _k \rho_k \mathcal N_x(\mu_k,\Sigma_k)$

Petersen y Pedersen dan la derivada parcial con respecto a la matriz de covarianza como$\Sigma$

$\frac{\delta \ln p(x)}{\delta \Sigma_j}=\frac{\rho_j\mathcal N_x(\mu_j,\Sigma_j)}{\sum _k \rho_k \mathcal N_x(\mu_k,\Sigma_k)}\frac{1}{2}[-\Sigma_j^{-1}+\Sigma_j^{-1}(x-\mu_j)(x-\mu_j)^T\Sigma_j^{-1}]$

El paso de descenso de gradiente para cada , como lo tengo implementado en Python es (esto es una ligera simplificación y el para todos los componentes se calcula antes de realizar la actualización):$\Sigma_j$$\Delta\Sigma$

j.sigma += learning_rate*(G(x)/M(x))*0.5*(-inv(j.sigma) + inv(j.sigma).dot((x-j.mu).dot((x-j.mu).transpose())).dot(inv(j.sigma)))

Donde j es un objeto que representa el componente de la mezcla y j.sigma y j.mu son la media y la varianza de ese componente. G (x) / M (x) a un código que calcula$j$$\frac{\rho_j\mathcal N_x(\mu_j,\Sigma_j)}{\sum _k \rho_k \mathcal N_x(\mu_k,\Sigma_k)}$

Entonces, me pregunto si hay algo mal con mi código (muy probable) o si esta es solo una forma realmente mala de ajustar este tipo de modelo cuando se trata de datos con más de dos dimensiones (consulte Toscano y McMurray para ver algoritmos para univariante y datos bivariados que definitivamente funcionan).

referencias: Toscano, JC y McMurray, B. (2010). Integración de señales con categorías: ponderación de señales acústicas en el habla utilizando aprendizaje no supervisado y estadísticas de distribución. Ciencia cognitiva, 34, 434-464.

Petersen y Pederson. The Matrix Cookbook, Versión: 14 de noviembre de 2008

Suponiendo que mus[d]es , es$\mu_j$j.sigma$\Sigma_j$, y de G(x)/M(x)hecho calcula la probabilidad posterior de componente$j$ dados los datos $x$,

• Esperaría que el acceso a la media, la covarianza y el cálculo de la parte posterior a todos impliquen cualquiera jo d, cualquier variable que represente el componente para el que desea calcular el gradiente en su código. Si nos dice qué jy qué drepresentamos, es posible que podamos brindarle más información.
• Si G(x)/M(x)accede j.Sigmaa calcular la parte posterior, es posible que su código no calcule lo que cree que hace. Sería mejor calcular primero todos los gradientes de todos los parámetros y luego realizar la actualización.

• El descenso de gradiente estocástico generalmente no es la primera opción para optimizar mezclas de gaussianos. Muy a menudo, se utiliza la maximización de expectativas (EM) (véase, por ejemplo, Bishop, 2007). Incluso si no usa EM, es posible que desee considerar BFGS o L-BFGS (implementado en scipy.optimize) antes de usar SGD. E incluso si se apega a SGD, debería considerar usar múltiples puntos de datos ("lotes") a la vez para estimar el gradiente, o al menos incluir un término de impulso . Mirando brevemente el artículo de Toscano y McMurray, supongo que eligieron usar SGD porque estaban interesados ​​en modelar la adquisición del habla de una manera biológicamente más plausible, en lugar de obtener el mejor ajuste posible, y hacerlo en línea (es decir, un dato punto a la vez). Si no necesita esto, mi consejo sería usar EM.

  (Me acabo de dar cuenta de que pediste específicamente el aprendizaje en línea , por lo que la única opción viable para ti podría ser agregar el término de impulso para acelerar un poco las cosas).

• La forma en que eligió calcular el gradiente es bastante ineficiente, lo que ralentizará aún más el aprendizaje. Es posible que no haya visto resultados razonables porque lleva mucho tiempo antes de que el algoritmo converja en algo interesante. Aquí hay una forma ligeramente mejor de calcular el gradiente:

  sigmaInv = inv(j.sigma)
  dSigma = G(x)/M(x) * 0.5 * (-sigmaInv + numpy.sum(sigmaInv.dot(x - mus[d]) * x))

  Todavía hay formas de mejorar aún más el cálculo del gradiente. Por ejemplo, todavía obtenemos una dirección de ascenso válida (aunque no un ascenso más pronunciado) si multiplicamos el gradiente por una matriz definida positiva (como$\Sigma_j$, lo que simplificaría un poco el gradiente). También podría funcionar mejor si utilizamos una parametrización diferente de la covarianza, como los factores de Cholesky , y calculamos los gradientes de esos.