Singularitätsprobleme im Gaußschen Mischungsmodell

In Kapitel 9 des Buches Mustererkennung und maschinelles Lernen gibt es diesen Teil über das Gaußsche Mischungsmodell:

Um ehrlich zu sein, verstehe ich nicht wirklich, warum dies eine Singularität schaffen würde. Kann mir das jemand erklären? Es tut mir leid, aber ich bin nur ein Student und ein Anfänger im maschinellen Lernen, daher mag meine Frage etwas albern klingen, aber bitte helfen Sie mir. Vielen Dank

Wenn wir einen Gaußschen Wert mit maximaler Wahrscheinlichkeit an einen einzelnen Datenpunkt anpassen möchten, erhalten wir einen sehr stacheligen Gaußschen Wert, der bis zu diesem Punkt "kollabiert". Die Varianz ist Null, wenn es nur einen Punkt gibt, was im Gaußschen Fall mit mehreren Variablen zu einer singulären Kovarianzmatrix führt, daher wird dies als Singularitätsproblem bezeichnet.

Wenn die Varianz Null wird, wird die Wahrscheinlichkeit der Gaußschen Komponente (Formel 9.15) unendlich und das Modell wird überpasst. Dies tritt nicht auf, wenn wir nur einen Gaußschen Wert an eine Anzahl von Punkten anpassen, da die Varianz nicht Null sein kann. Es kann jedoch vorkommen, dass wir eine Mischung aus Gaußschen haben, wie auf derselben Seite von PRML dargestellt.

Update :
Das Buch schlägt zwei Methoden zur Lösung des Singularitätsproblems vor:

1) Zurücksetzen des Mittelwerts und der Varianz, wenn Singularität auftritt

2) Verwenden von MAP anstelle von MLE durch Hinzufügen eines vorherigen.

Denken Sie daran, dass dieses Problem bei einer einzelnen Gaußschen Verteilung nicht aufgetreten ist. Um den Unterschied zu verstehen, ist zu beachten, dass, wenn ein einzelner Gaußscher auf einen Datenpunkt zusammenfällt, er multiplikative Faktoren zur Wahrscheinlichkeitsfunktion beiträgt, die sich aus den anderen Datenpunkten ergibt, und diese Faktoren exponentiell schnell auf Null gehen, was eine Gesamtwahrscheinlichkeit ergibt, die eher auf Null geht als die Unendlichkeit.

Ich bin auch ein bisschen durch diesen Teil verwirrt, und hier ist meine Interpretation. Nehmen Sie der Einfachheit halber 1D-Fall.

Wenn ein einzelner Gaußscher auf einem Datenpunkt , dh μ = x i , "kollabiert" , wird die Gesamtwahrscheinlichkeit:$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}} )$$

Sie sehen , wie , der Begriff auf der linken Seite p ( x i ) → ∞ , die wie die pathologischen Fall in GMM ist, aber der Begriff auf der rechten Seite , die die Wahrscheinlichkeit von anderen Datenpunkten ist p ( x ∖ i ) enthält noch Begriffe wie e - ( x n - μ ) $\sigma \to 0$$p(x_i) \to \infty$$p(\mathbf{x}\setminus{i})$ wobei→0exponentiell schnell ist wieσ→0, so dass der Gesamteffekt auf die Wahrscheinlichkeit darin besteht, dass es auf Null geht.$e^{-\frac{(x_n-\mu)^2}{2\sigma^2}}$$\to 0$$\sigma \to 0$

Der Hauptpunkt hierbei ist, dass beim Anpassen eines einzelnen Gaußschen Wertes alle Datenpunkte einen Parametersatz , anders als im Mischungsfall, in dem eine Komponente sich auf einen Datenpunkt "konzentrieren" kann, ohne die Gesamtdatenwahrscheinlichkeit zu beeinträchtigen .$\mu, \sigma$

Diese Antwort gibt einen Einblick in das Geschehen, das zu einer singulären Kovarianzmatrix beim Anpassen eines GMM an einen Datensatz führt, warum dies geschieht und was wir tun können, um dies zu verhindern.

Daher beginnen wir am besten damit, die Schritte während der Anpassung eines Gaußschen Mischungsmodells an einen Datensatz zusammenzufassen.

0. Entscheiden Sie, wie viele Quellen / Cluster (c) Sie an Ihre Daten
anpassen möchten. 1. Initialisieren Sie die Parameter Mittelwert , Kovarianz Σ c und Bruch_pro_Klasse π c pro Cluster c $\mu_c$$\Sigma_c$$\pi_c$

$\underline{E-Step}$

2. Berechnen Sie für jeden Datenpunkt die Wahrscheinlichkeit r i c, dass der Datenpunkt x i zum Cluster c gehört mit: r i c = π c N ( x i | μ c , Σ c$x_i$$r_{ic}$$x_i$
   
   wobeiN(x|μ,Σ)beschreibt den Gauß'schen mulitvariate mit: 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})$
   
   ricgibt uns für jeden Datenpunktxidas Maß:PRobeinbilitythatxibelongstoCla $$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$ Daherwennxiist sehr nahe an einem Gaußschen c, wird es ein hohes erhaltenricWert dafür gaußsche und sonst relativ niedrige Werte. M-Step_ Für jeden Cluster c: Berechnen Sie das Gesamtgewichtm$\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$(loses Sprechen des Bruchteils der Punkte, die dem Cluster c zugeordnet sind) und Aktualisieren von , μ c und Σ c unter Verwendung von r i c mit: 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}$$
   darandass Sie die aktualisierten Mittel in dieser letzten Formel verwenden. Iterativ wiederholen die E und MSchritt bis zur log-LikelihoodFunktion unseres Modell konvergiertwo das LogLikelihood mit berechnet wird: lnp(X|& pgr;,μ,Σ)=Σ N i = 1 ln(Σ $$\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}$$

ist nicht invertierbar und folgt Singular. Es ist auch plausibel, wenn wir davon ausgehen, dass die obige Matrix eine Matrix ist$A$ Es konnte keine Matrix geben $X$ was mit dieser Matrix punktiert die Identitätsmatrix gibt $I$(Nehmen Sie einfach diese Nullmatrix und addieren Sie sie mit jeder anderen 2x2-Matrix. Sie werden sehen, dass Sie immer die Nullmatrix erhalten.) Aber warum ist das ein Problem für uns? Betrachten Sie die obige Formel für die multivariate Norm. Dort würdest du finden$\boldsymbol{\Sigma_c^{-1}}$Das ist die Umkehrbarkeit der Kovarianzmatrix. Da eine singuläre Matrix nicht invertierbar ist, führt dies zu einem Fehler bei der Berechnung.
Nachdem wir nun wissen, wie eine singuläre, nicht invertierbare Matrix aussieht und warum dies für uns bei den GMM-Berechnungen wichtig ist, wie konnten wir auf dieses Problem stoßen? Zuallererst bekommen wir das$\boldsymbol{0}$Kovarianzmatrix oben, wenn der multivariate Gaußsche Wert während der Iteration zwischen dem E- und dem M-Schritt in einen Punkt fällt. Dies könnte passieren, wenn wir zum Beispiel einen Datensatz haben, zu dem wir 3 Gauß'sche Werte hinzufügen möchten, der aber eigentlich nur aus zwei Klassen (Clustern) besteht, so dass zwei dieser drei Gauß'schen Werte ihren eigenen Cluster abfangen, während der letzte Gauß'sche Wert nur diesen verwaltet einen einzigen Punkt zu fangen, auf dem es sitzt. Wir werden unten sehen, wie das aussieht. Aber Schritt für Schritt: Angenommen, Sie haben einen zweidimensionalen Datensatz, der aus zwei Clustern besteht, aber Sie wissen das nicht und möchten drei Gauß-Modelle anpassen, d. H. C = 3. Sie initialisieren Ihre Parameter im E-Schritt und zeichnen die Gaußschen oben auf Ihren Daten, die etw aussehen. wie (vielleicht können Sie die zwei relativ verstreuten Cluster links unten und rechts oben sehen): Nachdem Sie den Parameter initialisiert haben, führen Sie iterativ die Schritte E, T aus. Während dieses Vorgangs irren die drei Gaußschen herum und suchen nach ihrem optimalen Platz. Wenn Sie die Modellparameter beachten, ist das$\mu_c$ und $\pi_c$Sie werden feststellen, dass sie konvergieren, dass sie sich nach einer gewissen Anzahl von Iterationen nicht mehr ändern und damit der entsprechende Gaußsche seinen Platz im Raum gefunden hat. Wenn Sie eine Singularitätsmatrix haben, stoßen Sie auf etw. wie: Wo ich das dritte Gaußsche Modell mit rot eingekreist habe. Sie sehen also, dass dieser Gaußsche auf einem einzigen Datenpunkt sitzt, während die beiden anderen den Rest beanspruchen. Hier muss ich beachten, dass ich, um die Figur so zeichnen zu können, bereits die Kovarianz-Regularisierung verwendet habe, eine Methode zur Verhinderung von Singularitätsmatrizen, die im Folgenden beschrieben wird.

Ok, aber jetzt wissen wir immer noch nicht, warum und wie wir auf eine Singularitätsmatrix stoßen. Deshalb müssen wir uns die Berechnungen der$r_{ic}$ und die $cov$während der E- und M-Schritte. Wenn Sie sich das anschauen$r_{ic}$ formel nochmal:
$$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}})}$$ Sie sehen, dass dort die $r_{ic}$'s hätte große Werte, wenn sie sehr wahrscheinlich unter Cluster c sind, und niedrige Werte, wenn sie sonst sehr wahrscheinlich sind. Um dies deutlicher zu machen, betrachten wir den Fall, in dem wir zwei relativ verteilte Gaußsche und einen sehr engen Gaußschen haben und berechnen die$r_{ic}$ für jeden Datenpunkt $x_i$wie in der Abbildung dargestellt: Gehen Sie also die Datenpunkte von links nach rechts durch und stellen Sie sich vor, Sie würden die Wahrscheinlichkeit für jeden aufschreiben$x_i$dass es zum roten, blauen und gelben Gauß gehört. Was Sie sehen können, ist das für die meisten$x_i$Die Wahrscheinlichkeit, dass es zum gelben Gauß gehört, ist sehr gering. In dem obigen Fall, in dem der dritte Gaußsche auf einem einzelnen Datenpunkt sitzt,$r_{ic}$ ist nur für diesen einen Datenpunkt größer als Null, während es für jeden anderen Null ist $x_i$. (kollabiert auf diesen Datenpunkt -> Dies geschieht, wenn alle anderen Punkte wahrscheinlicher Teil von Gaußscher Eins oder Zwei sind und somit der einzige Punkt ist, der für Gaußsche Drei übrig bleibt -> Der Grund, warum dies geschieht, liegt in der Wechselwirkung zwischen der Datensatz selbst in der Initialisierung der Gaußschen. Wenn wir also andere Anfangswerte für die Gaußschen gewählt hätten, hätten wir ein anderes Bild gesehen, und die dritte Gaußsche würde möglicherweise nicht zusammenbrechen. Dies reicht aus, wenn Sie diesen Gaußschen immer weiter anstacheln. Das$r_{ic}$Tisch sieht dann aus wie: Wie Sie sehen können, die$r_{ic}$der dritten Spalte, das heißt für die dritte Gaußsche, sind Null anstelle dieser einen Zeile. Wenn wir nachsehen, welcher Datenpunkt hier dargestellt ist, erhalten wir den Datenpunkt: [23.38566343 8.07067598]. Ok, aber warum bekommen wir in diesem Fall eine Singularitätsmatrix? Nun, und dies ist unser letzter Schritt, deshalb müssen wir uns noch einmal mit der Berechnung der Kovarianzmatrix befassen, die lautet: $$\boldsymbol{\Sigma_c} \ = \ \Sigma_i r_{ic}(\boldsymbol{x_i}-\boldsymbol{\mu_c})^T(\boldsymbol{x_i}-\boldsymbol{\mu_c})$$ wir haben das alles gesehen $r_{ic}$ sind Null stattdessen für den einen $x_i$mit [23.38566343 8.07067598]. Jetzt soll die Formel von uns berechnet werden$(\boldsymbol{x_i}-\boldsymbol{\mu_c})$. Wenn wir uns das anschauen$\boldsymbol{\mu_c}$Für diesen dritten Gaußschen Wert erhalten wir [23.38566343 8.07067598]. Oh, aber warte, das ist genau das gleiche wie$x_i$ und so schrieb Bishop: "Nehmen wir an, eine der Komponenten des Mischungsmodells ist das $j$ Komponente hat ihren Mittelwert $\boldsymbol{\mu_j}$ genau gleich einem der Datenpunkte damit $\boldsymbol{\mu_j} = \boldsymbol{x_n}$für einen Wert von n "(Bishop, 2006, S.434). Also, was wird passieren? Nun, dieser Term wird Null sein, und daher war dieser Datenpunkt die einzige Chance für die Kovarianzmatrix, nicht Null zu werden (da dieser Datenpunkt Null war) der einzige wo$r_{ic}$> 0), es wird jetzt null und sieht so aus:

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

Folglich ist dies, wie oben erwähnt, eine singuläre Matrix und führt zu einem Fehler bei der Berechnung des multivariaten Gaußschen. Wie können wir eine solche Situation verhindern? Nun, wir haben gesehen, dass die Kovarianzmatrix singulär ist, wenn es die ist$\boldsymbol{0}$Matrix. Um Singularität zu verhindern, müssen wir einfach verhindern, dass die Kovarianzmatrix zu einer wird$\boldsymbol{0}$Matrix. Dazu wird der Digonalzahl der Kovarianzmatrix ein sehr kleiner Wert hinzugefügt (in sklearns GaussianMixture wird dieser Wert auf 1e-6 gesetzt). Es gibt auch andere Möglichkeiten, um Singularität zu verhindern, z. B. zu bemerken, wenn ein Gaußscher zusammenbricht, und seinen Mittelwert und / oder seine Kovarianzmatrix auf einen neuen, willkürlich hohen Wert einzustellen. Diese Kovarianz-Regularisierung ist auch im folgenden Code implementiert, mit dem Sie die beschriebenen Ergebnisse erhalten. Möglicherweise müssen Sie den Code mehrmals ausführen, um eine singuläre Kovarianzmatrix zu erhalten. Dies muss nicht jedes Mal geschehen, sondern hängt auch von der anfänglichen Einrichtung der Gaußschen ab.

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

Imho, alle Antworten verfehlen eine grundlegende Tatsache. Wenn man sich den Parameterraum für ein Gaußsches Mischungsmodell ansieht, ist dieser Raum entlang des Unterraums singulär, in dem weniger als die volle Anzahl von Komponenten in der Mischung vorhanden ist. Das bedeutet, dass Derivate automatisch Null sind und in der Regel der gesamte Unterraum als mle angezeigt wird. Aus philosophischer Sicht ist der Unterraum von weniger als vollständigen Kovarianzen die Grenze des Parameterraums, und man sollte immer misstrauisch sein, wenn das mle an der Grenze auftritt - es zeigt normalerweise an, dass ein größerer Parameterraum lauert, in dem man den finden kann "echte" mle. Es gibt ein Buch mit dem Titel "Algebraic Statistics" von Drton, Sturmfeld und Sullivant. Dieses Problem wird in diesem Buch ausführlich behandelt. Wenn Sie wirklich neugierig sind, sollten Sie sich das ansehen.

Für einen einzelnen Gaußschen Wert kann der Mittelwert möglicherweise einem der Datenpunkte entsprechen ($x_n$ zum Beispiel) und dann gibt es den folgenden Begriff in der Wahrscheinlichkeitsfunktion:

$$\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}$$ Das Limit $\sigma_j\rightarrow 0$ ist jetzt deutlich divergierend, da das Argument des Exponentials verschwindet.

Allerdings für einen Datenpunkt $x_m$ anders als der Mittelwert $\sigma_j$, wir werden haben

$$\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}$$ und jetzt divergiert das Argument des Exponentials (und ist negativ) in der Grenze $\sigma_j\rightarrow 0$. Infolgedessen verschwindet das Produkt dieser beiden Terme in der Wahrscheinlichkeitsfunktion.