Problemas de singularidad en el modelo de mezcla gaussiano

Question

En el capítulo 9 del libro Pattern recognition and machine learning, hay esta parte sobre el modelo de mezcla gaussiano:

Para ser honesto, no entiendo muy bien por qué esto crearía una singularidad. ¿Puede alguien explicármelo? Lo siento pero sólo soy un estudiante universitario y un novato en el aprendizaje automático, así que mi pregunta puede sonar un poco tonta, pero por favor ayúdenme. Muchas gracias

Answer 1

Si queremos ajustar una gaussiana a un único punto de datos utilizando la máxima verosimilitud, obtendremos una gaussiana muy puntiaguda que "colapsa" en ese punto. La varianza es cero cuando sólo hay un punto, lo que en el caso de la gaussiana multivariante conduce a una matriz de covarianza singular, por lo que se denomina problema de singularidad.

Cuando la varianza llega a cero, la probabilidad del componente gaussiano (fórmula 9.15) llega a infinito y el modelo se sobreajusta. Esto no ocurre cuando ajustamos sólo una gaussiana a un número de puntos, ya que la varianza no puede ser cero. Pero puede ocurrir cuando tenemos una mezcla de gaussianas, como se ilustra en la misma página de PRML.

Actualización :
El libro sugiere dos métodos para abordar el problema de la singularidad, que son

1) restablecer la media y la varianza cuando se produce una singularidad

2) utilizar MAP en lugar de MLE añadiendo una prioridad.

Answer 2

Esta respuesta le dará una idea de lo que está sucediendo que conduce a una matriz de covarianza singular durante el ajuste de un GMM a un conjunto de datos, por qué esto está sucediendo, así como lo que podemos hacer para evitar que.

Por lo tanto, lo mejor es empezar recapitulando los pasos durante el ajuste de un Modelo de Mezcla Gaussiana a un conjunto de datos.

0. Decide cuántas fuentes/clusters (c) quieres ajustar a tus datos
1. Inicializar los parámetros media $\mu_c$ covarianza $\Sigma_c$ y fracción_por_clase $\pi_c$ por clúster c

$\underline{E-Step}$

2. Calcular para cada punto de datos $x_i$ la probabilidad $r_{ic}$ ese dato $x_i$ pertenece al clúster c con:

   $$r_{ic} = \frac{\pi_c N(\boldsymbol{x_i} \ | \ \boldsymbol{\mu_c},\boldsymbol{\Sigma_c})}{\Sigma_{k=1}^K \pi_k N(\boldsymbol{x_i \ | \ \boldsymbol{\mu_k},\boldsymbol{\Sigma_k}})}$$ donde $N(\boldsymbol{x} \ | \ \boldsymbol{\mu},\boldsymbol{\Sigma})$ describe la gaussiana multivariante con:

   $$N(\boldsymbol{x_i},\boldsymbol{\mu_c},\boldsymbol{\Sigma_c}) \ = \ \frac{1}{(2\pi)^{\frac{n}{2}}|\Sigma_c|^{\frac{1}{2}}}exp(-\frac{1}{2}(\boldsymbol{x_i}-\boldsymbol{\mu_c})^T\boldsymbol{\Sigma_c^{-1}}(\boldsymbol{x_i}-\boldsymbol{\mu_c}))$$

   $r_{ic}$ nos da para cada punto de datos $x_i$ la medida de: $\frac{Probability \ that \ x_i \ belongs \ to \ class \ c}{Probability \ of \ x_i \ over \ all \ classes}$ por lo tanto, si $x_i$ está muy cerca de una gaussiana c, obtendrá un alto $r_{ic}$ para esta gaussiana y valores relativamente bajos en caso contrario.

   $\underline{M-Step}$

   Para cada conglomerado c: Calcular el peso total $m_c$ (en términos generales, la fracción de puntos asignados al conglomerado c) y actualizar $\pi_c$ , $\mu_c$ y $\Sigma_c$ utilizando $r_{ic}$ con:

   $$m_c \ = \ \Sigma_i r_ic$$
   $$\pi_c \ = \ \frac{m_c}{m}$$
   $$\boldsymbol{\mu_c} \ = \ \frac{1}{m_c}\Sigma_i r_{ic} \boldsymbol{x_i} $$
   $$\boldsymbol{\Sigma_c} \ = \ \frac{1}{m_c}\Sigma_i r_{ic}(\boldsymbol{x_i}-\boldsymbol{\mu_c})^T(\boldsymbol{x_i}-\boldsymbol{\mu_c})$$
   Ten en cuenta que en esta última fórmula tienes que utilizar los medios actualizados.

   Repite iterativamente los pasos E y M hasta que converja la función log-verosimilitud de nuestro modelo, donde la log-verosimilitud se calcula con:

   $$ln \ p(\boldsymbol{X} \ | \ \boldsymbol{\pi},\boldsymbol{\mu},\boldsymbol{\Sigma}) \ = \ \Sigma_{i=1}^N \ ln(\Sigma_{k=1}^K \pi_k N(\boldsymbol{x_i} \ | \ \boldsymbol{\mu_k},\boldsymbol{\Sigma_k}))$$

Así que ahora que hemos derivado los pasos singulares durante el cálculo tenemos que considerar lo que significa para una matriz ser singular. Una matriz es singular si no es invertible. Una matriz es invertible si existe una matriz $X$ tal que $AX = XA = I$ . Si no se da, se dice que la matriz es singular. Es decir, una matriz como:

\begin{bmatrix} 0 & 0 \\ 0 & 0 \end{bmatrix}

no es invertible y a continuación singular. También es plausible, que si suponemos que la matriz anterior es matriz $A$ no podría haber una matriz $X$ que da punteada con esta matriz la matriz identidad $I$ (Simplemente tome esta matriz cero y prodúzcala por puntos con cualquier otra matriz 2x2 y verá que siempre obtendrá la matriz cero). Pero, ¿por qué es esto un problema para nosotros? Bien, consideremos la fórmula para la normal multivariante anterior. En ella encontraríamos $\boldsymbol{\Sigma_c^{-1}}$ que es la invertible de la matriz de covarianza. Dado que una matriz singular no es invertible, esto nos arrojará un error durante el cálculo.
Ahora que sabemos cómo es una matriz singular no invertible y por qué es importante para nosotros durante los cálculos del MMG, ¿cómo podemos encontrarnos con este problema? En primer lugar, obtenemos lo siguiente $\boldsymbol{0}$ matriz de covarianza anterior si la gaussiana multivariante cae en un punto durante la iteración entre el paso E y M. Esto podría ocurrir si tenemos, por ejemplo, un conjunto de datos al que queremos ajustar 3 gaussianas, pero que en realidad sólo consta de dos clases (clusters), de modo que, hablando en términos generales, dos de estas tres gaussianas capturan su propio cluster, mientras que la última gaussiana sólo consigue capturar un único punto en el que se sitúa. Veremos cómo queda esto a continuación. Pero paso a paso: Supongamos que tenemos un conjunto de datos bidimensional que consta de dos conglomerados, pero no lo sabemos y queremos ajustarle tres modelos gaussianos, es decir, c = 3. Inicializamos los parámetros en el paso E y trazamos las gaussianas sobre los datos, que tienen un aspecto similar (tal vez podamos ver los dos conglomerados relativamente dispersos en la parte inferior izquierda y superior derecha): Una vez inicializado el parámetro, realiza iterativamente los pasos E, T. Durante este procedimiento, los tres gaussianos deambulan y buscan su lugar óptimo. Si se observan los parámetros del modelo, es decir $\mu_c$ y $\pi_c$ observarás que convergen, es decir, que después de un cierto número de iteraciones ya no cambian y con ello la gaussiana correspondiente ha encontrado su lugar en el espacio. En el caso de que usted tiene una matriz de singularidad se encuentra smth. como: Donde he rodeado de rojo el tercer modelo gaussiano. Para que veas que esta gaussiana se asienta en un único punto de datos, mientras que las otras dos reclaman el resto. Aquí tengo que notar que para ser capaz de dibujar la figura como que ya he utilizado la covarianza-regularización que es un método para prevenir las matrices de singularidad y se describe a continuación.

Ok , pero ahora todavía no sabemos por qué y cómo nos encontramos con una matriz de singularidad. Por lo tanto tenemos que mirar los cálculos de la $r_{ic}$ y el $cov$ durante los pasos E y M. Si se observa el $r_{ic}$ fórmula de nuevo:
$$r_{ic} = \frac{\pi_c N(\boldsymbol{x_i} \ | \ \boldsymbol{\mu_c},\boldsymbol{\Sigma_c})}{\Sigma_{k=1}^K \pi_k N(\boldsymbol{x_i \ | \ \boldsymbol{\mu_k},\boldsymbol{\Sigma_k}})}$$ ves que allí el $r_{ic}$ tendrían valores grandes si son muy probables bajo el conglomerado c y valores bajos en caso contrario. Para que esto resulte más evidente, consideremos el caso en el que tenemos dos gaussianas relativamente dispersas y una gaussiana muy ajustada y calculamos el $r_{ic}$ para cada punto de datos $x_i$ como se ilustra en la figura: Recorre los puntos de datos de izquierda a derecha e imagina que escribes la probabilidad de cada uno de ellos. $x_i$ que pertenece al gaussiano rojo, azul y amarillo. Lo que se ve es que para la mayoría de los $x_i$ la probabilidad de que pertenezca a la gaussiana amarilla es muy pequeña. En el caso anterior, en el que la tercera gaussiana se sitúa en un único punto de datos, $r_{ic}$ sólo es mayor que cero para este punto de datos, mientras que es cero para todos los demás. $x_i$ . (colapsa en este punto de datos --> Esto sucede si todos los demás puntos son más probablemente parte de la gaussiana uno o dos y por lo tanto este es el único punto que queda para la gaussiana tres --> La razón por la que esto sucede se puede encontrar en la interacción entre el conjunto de datos en sí en la inicialización de las gaussianas. Es decir, si hubiéramos elegido otros valores iniciales para las gaussianas, habríamos visto otra imagen y la tercera gaussiana quizás no colapsaría). Esto es suficiente si cada vez picas más esta gaussiana. El $r_{ic}$ mesa entonces parece smth. como: Como puede ver, el $r_{ic}$ de la tercera columna, es decir, para la tercera gaussiana son cero en lugar de esta fila. Si buscamos qué punto de datos está representado aquí obtenemos el punto de datos: [ 23.38566343 8.07067598]. Vale, pero ¿por qué obtenemos una matriz de singularidad en este caso? Bueno, y este es nuestro último paso, por lo tanto tenemos que considerar una vez más el cálculo de la matriz de covarianza que es: $$\boldsymbol{\Sigma_c} \ = \ \Sigma_i r_{ic}(\boldsymbol{x_i}-\boldsymbol{\mu_c})^T(\boldsymbol{x_i}-\boldsymbol{\mu_c})$$ hemos visto que todos $r_{ic}$ son cero en lugar de $x_i$ con [23,38566343 8,07067598]. Ahora la fórmula quiere que calculemos $(\boldsymbol{x_i}-\boldsymbol{\mu_c})$ . Si nos fijamos en el $\boldsymbol{\mu_c}$ para esta tercera gaussiana obtenemos [23,38566343 8,07067598]. Oh, pero espera, eso es exactamente lo mismo que $x_i$ y eso es lo que Bishop escribió con: "Supongamos que uno de los componentes del modelo de mezcla modelo, digamos el $j$ th componente, tiene su media $\boldsymbol{\mu_j}$ exactamente igual a uno de los puntos de datos para que $\boldsymbol{\mu_j} = \boldsymbol{x_n}$ para algún valor de n "(Bishop, 2006, p.434). ¿Qué ocurrirá entonces? Pues bien, este término será cero y, por lo tanto, este punto de datos era la única posibilidad de que la matriz de covarianza no fuera cero (ya que este punto de datos era el único en el que $r_{ic}$ >0), ahora se pone a cero y queda así:

\begin{bmatrix} 0 & 0 \\ 0 & 0 \end{bmatrix}

En consecuencia, como ya se ha dicho, se trata de una matriz singular y dará lugar a un error durante los cálculos de la gaussiana multivariante. Entonces, ¿cómo podemos evitar tal situación. Bien, hemos visto que la matriz de covarianza es singular si es la $\boldsymbol{0}$ matriz. Por lo tanto para evitar la singularidad simplemente tenemos que evitar que la matriz de covarianza se convierta en una $\boldsymbol{0}$ matriz. Esto se hace añadiendo un valor muy pequeño (en GaussianMixture de sklearn este valor se fija en 1e-6) a la digonal de la matriz de covarianza. También hay otras formas de evitar la singularidad, como darse cuenta de cuándo colapsa una gaussiana y fijar su media y/o su matriz de covarianza en un nuevo valor o valores arbitrariamente altos. Esta regularización de covarianza también se implementa en el código de abajo con el que se obtienen los resultados descritos. Tal vez usted tiene que ejecutar el código varias veces para obtener una matriz de covarianza singular, ya que, como se ha dicho. esto no debe suceder cada vez, sino que también depende de la configuración inicial de los gaussianos.

import matplotlib.pyplot as plt
from matplotlib import style
style.use('fivethirtyeight')
from sklearn.datasets.samples_generator import make_blobs
import numpy as np
from scipy.stats import multivariate_normal

# 0. Create dataset
X,Y = make_blobs(cluster_std=2.5,random_state=20,n_samples=500,centers=3)

# Stratch dataset to get ellipsoid data
X = np.dot(X,np.random.RandomState(0).randn(2,2))

class EMM:

    def __init__(self,X,number_of_sources,iterations):
        self.iterations = iterations
        self.number_of_sources = number_of_sources
        self.X = X
        self.mu = None
        self.pi = None
        self.cov = None
        self.XY = None

    # Define a function which runs for i iterations:
    def run(self):
        self.reg_cov = 1e-6*np.identity(len(self.X[0]))
        x,y = np.meshgrid(np.sort(self.X[:,0]),np.sort(self.X[:,1]))
        self.XY = np.array([x.flatten(),y.flatten()]).T

        # 1. Set the initial mu, covariance and pi values
        self.mu = np.random.randint(min(self.X[:,0]),max(self.X[:,0]),size=(self.number_of_sources,len(self.X[0]))) # This is a nxm matrix since we assume n sources (n Gaussians) where each has m dimensions
        self.cov = np.zeros((self.number_of_sources,len(X[0]),len(X[0]))) # We need a nxmxm covariance matrix for each source since we have m features --> We create symmetric covariance matrices with ones on the digonal
        for dim in range(len(self.cov)):
            np.fill_diagonal(self.cov[dim],5)

        self.pi = np.ones(self.number_of_sources)/self.number_of_sources # Are "Fractions"
        log_likelihoods = [] # In this list we store the log likehoods per iteration and plot them in the end to check if
                             # if we have converged

        # Plot the initial state    
        fig = plt.figure(figsize=(10,10))
        ax0 = fig.add_subplot(111)
        ax0.scatter(self.X[:,0],self.X[:,1])
        for m,c in zip(self.mu,self.cov):
            c += self.reg_cov
            multi_normal = multivariate_normal(mean=m,cov=c)
            ax0.contour(np.sort(self.X[:,0]),np.sort(self.X[:,1]),multi_normal.pdf(self.XY).reshape(len(self.X),len(self.X)),colors='black',alpha=0.3)
            ax0.scatter(m[0],m[1],c='grey',zorder=10,s=100)

        mu = []
        cov = []
        R = []

        for i in range(self.iterations):               

            mu.append(self.mu)
            cov.append(self.cov)

            # E Step
            r_ic = np.zeros((len(self.X),len(self.cov)))

            for m,co,p,r in zip(self.mu,self.cov,self.pi,range(len(r_ic[0]))):
                co+=self.reg_cov
                mn = multivariate_normal(mean=m,cov=co)
                r_ic[:,r] = p*mn.pdf(self.X)/np.sum([pi_c*multivariate_normal(mean=mu_c,cov=cov_c).pdf(X) for pi_c,mu_c,cov_c in zip(self.pi,self.mu,self.cov+self.reg_cov)],axis=0)
            R.append(r_ic)

            # M Step

            # Calculate the new mean vector and new covariance matrices, based on the probable membership of the single x_i to classes c --> r_ic
            self.mu = []
            self.cov = []
            self.pi = []
            log_likelihood = []

            for c in range(len(r_ic[0])):
                m_c = np.sum(r_ic[:,c],axis=0)
                mu_c = (1/m_c)*np.sum(self.X*r_ic[:,c].reshape(len(self.X),1),axis=0)
                self.mu.append(mu_c)

                # Calculate the covariance matrix per source based on the new mean
                self.cov.append(((1/m_c)*np.dot((np.array(r_ic[:,c]).reshape(len(self.X),1)*(self.X-mu_c)).T,(self.X-mu_c)))+self.reg_cov)
                # Calculate pi_new which is the "fraction of points" respectively the fraction of the probability assigned to each source 
                self.pi.append(m_c/np.sum(r_ic)) 

            # Log likelihood
            log_likelihoods.append(np.log(np.sum([k*multivariate_normal(self.mu[i],self.cov[j]).pdf(X) for k,i,j in zip(self.pi,range(len(self.mu)),range(len(self.cov)))])))

        fig2 = plt.figure(figsize=(10,10))
        ax1 = fig2.add_subplot(111) 
        ax1.plot(range(0,self.iterations,1),log_likelihoods)
        #plt.show()
        print(mu[-1])
        print(cov[-1])
        for r in np.array(R[-1]):
            print(r)
        print(X)

    def predict(self):
        # PLot the point onto the fittet gaussians
        fig3 = plt.figure(figsize=(10,10))
        ax2 = fig3.add_subplot(111)
        ax2.scatter(self.X[:,0],self.X[:,1])
        for m,c in zip(self.mu,self.cov):
            multi_normal = multivariate_normal(mean=m,cov=c)
            ax2.contour(np.sort(self.X[:,0]),np.sort(self.X[:,1]),multi_normal.pdf(self.XY).reshape(len(self.X),len(self.X)),colors='black',alpha=0.3)

EMM = EMM(X,3,100)     
EMM.run()
EMM.predict()

Answer 3

Recordemos que este problema no se planteó en el distribución gaussiana. Para entender la diferencia, obsérvese que si a gaussiana colapsa en un punto de datos, contribuirá con factores factores multiplicativos a la función de verosimilitud derivada de las otros puntos de datos y estos factores llegarán a cero exponencialmente exponencialmente rápido, dando una probabilidad global que va a cero en lugar de al infinito. infinito.

Yo también estoy algo confuso con esta parte, y ésta es mi interpretación. Tomemos el caso 1D para simplificar.

Cuando una sola gaussiana "colapsa" sobre un punto de datos $x_i$ es decir, $\mu=x_i$ la probabilidad total se convierte en:

$$p(\mathbf{x}) = p(x_i) p(\mathbf{x}\setminus{i}) = (\frac{1}{\sqrt{2\pi}\sigma}) (\prod_{n \neq i}^N \frac{1}{\sqrt{2\pi}\sigma} e^{-\frac{(x_n-\mu)^2}{2\sigma^2}} )$$

Usted ve como $\sigma \to 0$ el término de la izquierda $p(x_i) \to \infty$ que es como el caso patológico en GMM, pero el término de la derecha, que es la probabilidad de otros puntos de datos $p(\mathbf{x}\setminus{i})$ todavía contiene términos como $e^{-\frac{(x_n-\mu)^2}{2\sigma^2}}$ que $\to 0$ exponencialmente rápido como $\sigma \to 0$ por lo que el efecto global sobre la probabilidad es que vaya a cero.

El punto principal aquí es que al ajustar una sola gaussiana, todos los puntos de datos tienen que compartir un conjunto de parámetros $\mu, \sigma$ a diferencia del caso de las mezclas, en el que un componente puede "centrarse" en un punto de datos sin afectar a la probabilidad global de los datos.

Answer 4

En mi opinión, todas las respuestas pasan por alto un hecho fundamental. Si se observa el espacio de parámetros de un modelo de mezcla gaussiano, este espacio es singular a lo largo del subespacio en el que hay menos del número total de componentes de la mezcla. Esto significa que las derivadas son automáticamente cero y normalmente todo el subespacio se mostrará como un mle. Más filosóficamente, el subespacio de menos de covarianzas de rango completo es el límite del espacio de parámetros y uno siempre debe sospechar cuando el mle se produce en el límite - por lo general indica que hay un mayor espacio de parámetros al acecho en el que uno puede encontrar el "real" mle. Hay un libro llamado "Algebraic Statistics" por Drton, Sturmfeld, y Sullivant. En ese libro se trata esta cuestión con cierto detalle. Si usted es realmente curioso, usted debe mirar eso.

Answer 5

Para una gaussiana simple, la media puede ser igual a uno de los puntos de datos ( $x_n$ por ejemplo) y entonces aparece el siguiente término en la función de verosimilitud:

\begin{equation} {\cal N}(x_n|x_n,\sigma_j 1\!\!1)\rightarrow \lim_{\sigma_j\rightarrow x_n}\frac{1}{(2\pi)^{1/2}\sigma_j} \exp \left( -\frac{1}{\sigma_j}|x_n-\sigma_j|^2 \right)= \frac{1}{(2\pi)^{1/2}\sigma_j} \end{equation} El límite $\sigma_j\rightarrow 0$ es ahora claramente divergente ya que el argumento de la exponencial desaparece.

Sin embargo, para un punto de datos $x_m$ diferente de la media $\sigma_j$ tendremos \begin{equation} {\cal N}(x_m|x_m,\sigma_j 1\!\!1)= \frac{1}{(2\pi)^{1/2}\sigma_j} \exp \left( -\frac{1}{\sigma_j}|x_m-\sigma_j|^2 \right) \end{equation} y ahora el argumento de la exponencial diverge (y es negativo) en el límite $\sigma_j\rightarrow 0$ . Como resultado, el producto de estos dos términos en la función de probabilidad desaparecerá.