Проблемы сингулярности в модели гауссовой смеси

В главе 9 книги «Распознавание образов и машинное обучение» описана модель гауссовой смеси:

Честно говоря, я не очень понимаю, почему это создаст особенность. Кто-нибудь может мне это объяснить? Извините, но я всего лишь студент и новичок в машинном обучении, поэтому мой вопрос может показаться немного глупым, но, пожалуйста, помогите мне. большое спасибо

Если мы хотим подогнать гауссиан к одной точке данных, используя максимальную вероятность, мы получим очень колючий гауссиан, который «сжимается» до этой точки. Дисперсия равна нулю, когда есть только одна точка, которая в многовариантном гауссовском случае приводит к сингулярной ковариационной матрице, поэтому она называется проблемой сингулярности.

Когда дисперсия достигает нуля, вероятность гауссовой составляющей (формула 9.15) уходит в бесконечность, и модель становится переобученной. Этого не происходит, когда мы подгоняем только один гауссов к числу точек, поскольку дисперсия не может быть нулевой. Но это может произойти, когда у нас есть смесь гауссов, как показано на той же странице PRML.

Обновление :
книга предлагает два метода решения проблемы сингулярности, которые

1) сброс среднего значения и дисперсии при возникновении сингулярности

2) используя MAP вместо MLE, добавив априор.

Напомним, что эта проблема не возникла в случае одного гауссовского распределения. Чтобы понять разницу, обратите внимание, что если один гауссовский коллапс падает на точку данных, он будет вносить мультипликативные факторы в функцию правдоподобия, возникающую из других точек данных, и эти коэффициенты будут стремительно стремиться к нулю в геометрической прогрессии, что дает общую вероятность, которая скорее достигает нуля чем бесконечность.

Я также немного смущен этой частью, и вот моя интерпретация. Возьмите 1D случай для простоты.

Когда один гауссов «коллапсирует» в точке данных , то есть μ = x i , общая вероятность становится равной:$x_i$$\mu=x_i$

$$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}} )$$

Вы видите , как , то член в левом р ( х я ) → ∞ , которая, как патологический случай в GMM, но член в правой части , которая является вероятностью других точек данных р ( х ∖ я ) , по-прежнему содержит такие термины, как e - ( x n - μ ) $\sigma \to 0$$p(x_i) \to \infty$$p(\mathbf{x}\setminus{i})$ который→0экспоненциально быстро приσ→0, поэтому общее влияние на вероятность состоит в том, что оно стремится к нулю.$e^{-\frac{(x_n-\mu)^2}{2\sigma^2}}$$\to 0$$\sigma \to 0$

Основной момент здесь заключается в том, что при подборе одного гауссиана все точки данных должны совместно использовать один набор параметров , в отличие от случая смешения, когда один компонент может «фокусироваться» на одной точке данных без ущерба для общей вероятности данных ,$\mu, \sigma$

Этот ответ даст представление о том, что происходит, что приводит к особой ковариационной матрице во время подгонки GMM к набору данных, почему это происходит, а также о том, что мы можем сделать, чтобы предотвратить это.

Поэтому лучше всего начать с повторения этапов во время подгонки модели гауссовой смеси к набору данных.

0. Определите, сколько источников / кластеров (c) вы хотите разместить в своих данных.
1. Инициализируйте средние значения параметров , ковариацию Σ c и Fraction_per_class π c на кластер c. $\mu_c$$\Sigma_c$$\pi_c$

$\underline{E-Step}$

2. Рассчитать для каждой точки данных вероятностной г я гр что DataPoint х я принадлежу кластеру с с: г я с = π с N ( х я | ц с , Σ с$x_i$$r_{ic}$$x_i$
   
   гдеN(x|μ,Σ)описывает многовариантный гауссов с: N(xi,μc,Σc)= $$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}})}$$ $N(\boldsymbol{x} \ | \ \boldsymbol{\mu},\boldsymbol{\Sigma})$
   
   ricдает нам для каждого элемента данныхxiмеру:Probabilitythatхябелонгстосл $$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}$$x_i$ , следовательноеслихяочень близко к одному гауссовой с, он будет получать высокуюRIсзначением для этого гауссова и относительно низкие значения в противном случае. M-Step_ Для каждого кластера c: Рассчитать общий весm$\frac{Probability \ that \ x_i \ belongs \ to \ class \ c}{Probability \ of \ x_i \ over \ all \ classes}$$x_i$$r_{ic}$
   
   $\underline{M-Step}$
   
   $m_c$(условно говоря, доля точек, выделенных кластеру c) и обновите , µ c и Σ c, используя r i c, с помощью: m c = Σ i r i c π c = m $\pi_c$$\mu_c$$\Sigma_c$$r_{ic}$
   
   $$m_c \ = \ \Sigma_i r_ic$$
   μc= $$\pi_c \ = \ \frac{m_c}{m}$$
   Σc= $$\boldsymbol{\mu_c} \ = \ \frac{1}{m_c}\Sigma_i r_{ic} \boldsymbol{x_i}$$
   Помните, что вы должны использовать обновленные средства в этой последней формуле. Итеративно повторяют Е и М шагпока функции логарифмического правдоподобия нашей модели сходитсягде журнал правдоподобия вычисляется с: лпр(Х|П,ц,Σ)=Е N я = 1 лп(Е $$\boldsymbol{\Sigma_c} \ = \ \frac{1}{m_c}\Sigma_i r_{ic}(\boldsymbol{x_i}-\boldsymbol{\mu_c})^T(\boldsymbol{x_i}-\boldsymbol{\mu_c})$$
   
   
   
   
   $$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}))$$

$X$$AX = XA = I$

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

$A$$X$$I$$\boldsymbol{\Sigma_c^{-1}}$
$\boldsymbol{0}$ковариационная матрица выше, если многовариантный гауссов попадает в одну точку во время итерации между этапами E и M. Это может произойти, если у нас есть, например, набор данных, к которому мы хотим вписать 3 гауссиана, но который на самом деле состоит только из двух классов (кластеров), так что, грубо говоря, два из этих трех гауссианов ловят свой кластер, а последний гауссиан управляет им только поймать одну единственную точку, на которой он сидит. Посмотрим, как это будет выглядеть ниже. Но шаг за шагом: предположим, что у вас есть двумерный набор данных, который состоит из двух кластеров, но вы этого не знаете и хотите подогнать к нему три гауссовых модели, то есть c = 3. Вы инициализируете свои параметры в шаге E и строите график гауссиане на вершине ваших данных, которые выглядят чем-то. например (возможно, вы можете видеть два относительно рассеянных кластера слева внизу и вверху справа): $\mu_c$$\pi_c$

$r_{ic}$$cov$$r_{ic}$
$$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}})}$$ $r_{ic}$$r_{ic}$$x_i$$x_i$$x_i$$r_{ic}$$x_i$$r_{ic}$$r_{ic}$ $$\boldsymbol{\Sigma_c} \ = \ \Sigma_i r_{ic}(\boldsymbol{x_i}-\boldsymbol{\mu_c})^T(\boldsymbol{x_i}-\boldsymbol{\mu_c})$$ $r_{ic}$$x_i$$(\boldsymbol{x_i}-\boldsymbol{\mu_c})$$\boldsymbol{\mu_c}$$x_i$$j$$\boldsymbol{\mu_j}$$\boldsymbol{\mu_j} = \boldsymbol{x_n}$$r_{ic}$

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

$\boldsymbol{0}$$\boldsymbol{0}$матрица. Это делается путем добавления очень небольшого значения (в GaussianMixture у sklearn это значение равно 1e-6) к дигональному ковариационной матрицы. Существуют также другие способы предотвращения сингулярности, такие как уведомление о коллапсе гауссиана и установка его среднего значения и / или ковариационной матрицы на новое, произвольно высокое значение (я). Эта ковариационная регуляризация также реализована в приведенном ниже коде, с помощью которого вы получите описанные результаты. Возможно, вам придется запускать код несколько раз, чтобы получить единственную ковариационную матрицу, как уже было сказано. это не должно происходить каждый раз, но также зависит от первоначальной настройки гауссиан.

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

Имхо, во всех ответах пропущен фундаментальный факт. Если взглянуть на пространство параметров для модели гауссовой смеси, это пространство является сингулярным вдоль подпространства, где в смеси меньше полного числа компонентов. Это означает, что производные автоматически равны нулю, и обычно все подпространство будет отображаться как mle. Говоря более философски, подпространство ковариаций менее полного ранга является границей пространства параметров, и всегда следует быть подозрительным, когда на границе возникает mle - обычно это указывает на то, что вокруг скрывается большее пространство параметров, в котором можно найти "настоящий" мле. Есть книга под названием «Алгебраическая статистика» Дртона, Штурмфельда и Салливанта. Эта проблема обсуждается в этой книге более подробно. Если вам действительно любопытно, вы должны посмотреть на это.

$x_n$

$$\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}$$ Предел $\sigma_j\rightarrow 0$ теперь явно расходится, так как аргумент экспоненты исчезает.

Однако для точки данных $x_m$ отличается от среднего $\sigma_j$, у нас будет

$$\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}$$ и теперь аргумент экспоненты расходится (и является отрицательным) в пределе $\sigma_j\rightarrow 0$, В результате произведение этих двух членов в функции правдоподобия исчезнет.