9.1 The Gaussian mixture model and KK-means

We return once again to the GMM, but this time, armed with the EM algorithm, we finally derive the learning rules. The joint cross entropy for a Gaussian mixture model is

H(p⁒pΛ‡)⁒p^⁒[𝑿ˇ,𝒀;𝜽]β‰ˆβŸ¨βˆ’log⁑p^⁒(𝑿ˇ,𝒀;𝜽)βŸ©π‘ΏΛ‡,𝒀=βŸ¨βˆ’logp^(𝒀|𝑿ˇ;𝜽)p^(𝑿ˇ;𝜽)βŸ©π‘ΏΛ‡,𝒀=βŸ¨βˆ’log⁒∏k=1K[𝒩︀⁒(𝒀;𝝁k,𝚺k)⁒πk]XΛ‡kβŸ©π‘ΏΛ‡,𝒀=+βŸ¨βˆ’βˆ‘k=1KXΛ‡k⁒[12⁒log⁑|𝚺kβˆ’1|βˆ’12⁒(𝝁kβˆ’π’€)T⁒𝚺kβˆ’1⁒(𝝁kβˆ’π’€)+log⁑πk]βŸ©π‘ΏΛ‡,𝒀.\begin{split}{\text{H}_{(p\check{p})\hat{p}}{\mathopen{}\mathclose{{}\left[{% \bm{\check{X}}}{},{\bm{Y}}{};\bm{\theta}}\right]}}&{}\approx{\mathopen{}% \mathclose{{}\left\langle{-\log{\hat{p}\mathopen{}\mathclose{{}\left({\bm{% \check{X}}},{\bm{Y}};\bm{\theta}}\right)}}}\right\rangle_{{\bm{\check{X}}}{},{% \bm{Y}}{}}}\\ &{}={\mathopen{}\mathclose{{}\left\langle{-\log{\hat{p}\mathopen{}\mathclose{{% }\left({\bm{Y}}\middle|{\bm{\check{X}}};\bm{\theta}}\right)}{\hat{p}\mathopen{% }\mathclose{{}\left({\bm{\check{X}}};\bm{\theta}}\right)}}}\right\rangle_{{\bm% {\check{X}}}{},{\bm{Y}}{}}}\\ &{}={\mathopen{}\mathclose{{}\left\langle{-\log\prod_{k=1}^{{K}}\mathopen{}% \mathclose{{}\left[\mathcal{N}\mathopen{}\mathclose{{}\left({\bm{Y}};\bm{\mu}_% {k},\>\mathbf{{\Sigma}}_{k}}\right)\pi_{k}}\right]^{{\check{X}}_{k}}}}\right% \rangle_{{\bm{\check{X}}}{},{\bm{Y}}{}}}\\ &{}\stackrel{{\scriptstyle\text{+}}}{{=}}{\mathopen{}\mathclose{{}\left\langle% {-\sum_{k=1}^{{K}}{\check{X}}_{k}\mathopen{}\mathclose{{}\left[\frac{1}{2}\log% \mathopen{}\mathclose{{}\left\lvert\mathbf{{\Sigma}}_{k}^{-1}}\right\rvert-% \frac{1}{2}\mathopen{}\mathclose{{}\left(\bm{\mu}_{k}-{\bm{Y}}}\right)^{\text{% T}}\mathbf{{\Sigma}}_{k}^{-1}\mathopen{}\mathclose{{}\left(\bm{\mu}_{k}-{\bm{Y% }}}\right)+\log\pi_{k}}\right]}}\right\rangle_{{\bm{\check{X}}}{},{\bm{Y}}{}}}% .\end{split}

To enforce the fact that the prior probabilities sum to one, we can augment the loss with a Lagrangian term:

ℒ︀⁒(𝜽)=λ⁒(βˆ‘k=1KΟ€kβˆ’1)+H(p⁒pΛ‡)⁒p^⁒[𝑿ˇ,𝒀;𝜽]=+λ⁒(βˆ‘k=1KΟ€kβˆ’1)+βŸ¨βˆ’βˆ‘k=1KXΛ‡k⁒[12⁒log⁑|𝚺kβˆ’1|βˆ’12⁒(𝝁kβˆ’π’€)T⁒𝚺kβˆ’1⁒(𝝁kβˆ’π’€)+log⁑πk]βŸ©π‘ΏΛ‡,𝒀.\begin{split}\mathcal{L}(\bm{\theta})&{}=\lambda\mathopen{}\mathclose{{}\left(% \sum_{k=1}^{{K}}\pi_{k}-1}\right)+{\text{H}_{(p\check{p})\hat{p}}{\mathopen{}% \mathclose{{}\left[{\bm{\check{X}}}{},{\bm{Y}}{};\bm{\theta}}\right]}}\\ &{}\stackrel{{\scriptstyle\text{+}}}{{=}}\lambda\mathopen{}\mathclose{{}\left(% \sum_{k=1}^{{K}}\pi_{k}-1}\right)+{\mathopen{}\mathclose{{}\left\langle{-\sum_% {k=1}^{{K}}{\check{X}}_{k}\mathopen{}\mathclose{{}\left[\frac{1}{2}\log% \mathopen{}\mathclose{{}\left\lvert\mathbf{{\Sigma}}_{k}^{-1}}\right\rvert-% \frac{1}{2}\mathopen{}\mathclose{{}\left(\bm{\mu}_{k}-{\bm{Y}}}\right)^{\text{% T}}\mathbf{{\Sigma}}_{k}^{-1}\mathopen{}\mathclose{{}\left(\bm{\mu}_{k}-{\bm{Y% }}}\right)+\log\pi_{k}}\right]}}\right\rangle_{{\bm{\check{X}}}{},{\bm{Y}}{}}}% .\end{split}

The M step.

We take the derivatives in turn. First the mixing proportions:

0=setd⁒ℒ︀d⁒πk=Ξ»βˆ’βŸ¨XΛ‡kβŸ©π‘ΏΛ‡,𝒀πkβŸΉβˆ‘k=1KΟ€k⁒λ=βˆ‘k=1K⟨XΛ‡kβŸ©π‘ΏΛ‡,π’€βŸΉΞ»=1βŸΉΟ€k=⟨XΛ‡kβŸ©π‘ΏΛ‡,𝒀;0\stackrel{{\scriptstyle\text{set}}}{{=}}\frac{\mathrm{d}{\mathcal{L}}}{% \mathrm{d}{\pi_{k}}}=\lambda-\frac{{\mathopen{}\mathclose{{}\left\langle{{% \check{X}}_{k}}}\right\rangle_{{\bm{\check{X}}}{},{\bm{Y}}{}}}}{\pi_{k}}% \implies\sum_{k=1}^{{K}}\pi_{k}\lambda=\sum_{k=1}^{{K}}{\mathopen{}\mathclose{% {}\left\langle{{\check{X}}_{k}}}\right\rangle_{{\bm{\check{X}}}{},{\bm{Y}}{}}}% \implies\lambda=1\implies\pi_{k}={\mathopen{}\mathclose{{}\left\langle{{\check% {X}}_{k}}}\right\rangle_{{\bm{\check{X}}}{},{\bm{Y}}{}}};

then the class-conditional means:

0=setd⁒ℒ︀d⁒𝝁kT=⟨XΛ‡k⁒(𝝁kβˆ’π’€)T⁒𝚺kβˆ’1βŸ©π‘ΏΛ‡,π’€βŸΉπk=⟨XΛ‡kβ’π’€βŸ©π‘ΏΛ‡,π’€βŸ¨XΛ‡kβŸ©π‘ΏΛ‡,𝒀;0\stackrel{{\scriptstyle\text{set}}}{{=}}\frac{\mathrm{d}{\mathcal{L}}}{% \mathrm{d}{\bm{\mu}_{k}}^{\text{T}}}={\mathopen{}\mathclose{{}\left\langle{{% \check{X}}_{k}\mathopen{}\mathclose{{}\left(\bm{\mu}_{k}-{\bm{Y}}}\right)^{% \text{T}}\mathbf{{\Sigma}}_{k}^{-1}}}\right\rangle_{{\bm{\check{X}}}{},{\bm{Y}% }{}}}\implies\bm{\mu}_{k}=\frac{{\mathopen{}\mathclose{{}\left\langle{{\check{% X}}_{k}{\bm{Y}}}}\right\rangle_{{\bm{\check{X}}}{},{\bm{Y}}{}}}}{{\mathopen{}% \mathclose{{}\left\langle{{\check{X}}_{k}}}\right\rangle_{{\bm{\check{X}}}{},{% \bm{Y}}{}}}};

and the class-conditional covariances:

0=setd⁒ℒ︀d⁒𝚺kβˆ’1=βŸ¨βˆ’XΛ‡k⁒[12⁒𝚺kβˆ’12⁒(𝝁kβˆ’π’€)⁒(𝝁kβˆ’π’€)T]βŸ©π‘ΏΛ‡,π’€βŸΉπšΊk=⟨XΛ‡k⁒(𝝁kβˆ’π’€)⁒(𝝁kβˆ’π’€)TβŸ©π‘ΏΛ‡,π’€βŸ¨XΛ‡kβŸ©π‘ΏΛ‡,𝒀=⟨XΛ‡k⁒𝒀⁒𝒀TβŸ©π‘ΏΛ‡,π’€βŸ¨XΛ‡kβŸ©π‘ΏΛ‡,π’€βˆ’πk⁒𝝁kT.\begin{split}0\stackrel{{\scriptstyle\text{set}}}{{=}}\frac{\mathrm{d}{% \mathcal{L}}}{\mathrm{d}{\mathbf{{\Sigma}}_{k}^{-1}}}&{}={\mathopen{}% \mathclose{{}\left\langle{-{\check{X}}_{k}\mathopen{}\mathclose{{}\left[\frac{% 1}{2}\mathbf{{\Sigma}}_{k}-\frac{1}{2}\mathopen{}\mathclose{{}\left(\bm{\mu}_{% k}-{\bm{Y}}}\right)\mathopen{}\mathclose{{}\left(\bm{\mu}_{k}-{\bm{Y}}}\right)% ^{\text{T}}}\right]}}\right\rangle_{{\bm{\check{X}}}{},{\bm{Y}}{}}}\\ \implies\mathbf{{\Sigma}}_{k}&{}=\frac{{\mathopen{}\mathclose{{}\left\langle{{% \check{X}}_{k}\mathopen{}\mathclose{{}\left(\bm{\mu}_{k}-{\bm{Y}}}\right)% \mathopen{}\mathclose{{}\left(\bm{\mu}_{k}-{\bm{Y}}}\right)^{\text{T}}}}\right% \rangle_{{\bm{\check{X}}}{},{\bm{Y}}{}}}}{{\mathopen{}\mathclose{{}\left% \langle{{\check{X}}_{k}}}\right\rangle_{{\bm{\check{X}}}{},{\bm{Y}}{}}}}=\frac% {{\mathopen{}\mathclose{{}\left\langle{{\check{X}}_{k}{\bm{Y}}{\bm{Y}}^{\text{% T}}}}\right\rangle_{{\bm{\check{X}}}{},{\bm{Y}}{}}}}{{\mathopen{}\mathclose{{}% \left\langle{{\check{X}}_{k}}}\right\rangle_{{\bm{\check{X}}}{},{\bm{Y}}{}}}}-% \bm{\mu}_{k}\bm{\mu}_{k}^{\text{T}}.\end{split}

Whether in EM or under a fully observed model, the optimal parameters are intuitive. The optimal mixing proportion, emission mean, and emission covariance for class kk are their sample counterparts, i.e.Β the sample proportion, sample mean, and sample covariance (resp.); or, to put it yet another way, the average number of occurrences of class kk, the average value of the samples 𝒀{\bm{Y}} from class kk, and the average covariance of the samples 𝒀{\bm{Y}} from class kk. The difference between learning (in EM) with latent, rather than fully observed, classes is that these are weighted, rather than unweighted, averages. In particular, class assignments are soft:Β each class takes some continuous-valued responsibility††margin: class responsibilities for each datum π’šn\bm{y}_{n}, namely 𝔼p^[XΛ‡k|π’šn]=p^(XΛ‡k=1|π’šn;𝜽old)\mathbb{E}_{\hat{p}}{\mathopen{}\mathclose{{}\left[{\check{X}}_{k}|\bm{y}_{n}}% \right]}={\hat{p}\mathopen{}\mathclose{{}\left({\check{X}}_{k}=1\middle|\bm{y}% _{n};\bm{\theta}_{\text{old}}}\right)}, the probability of that class under the (previous) posterior distribution.

For example, when the class labels are observed, the denominator in the equation for the optimal mean becomes just the number of times class kk occurred, and the numerator picks out just those samples π’šn\bm{y}_{n} associated with class kk. This is the sample average of the π’šn\bm{y}_{n} in class kk. But when the class 𝑿ˇ{\bm{\check{X}}} is latent, the denominator is the average soft assignment or responsibility of class kk, and the numerator is the (soft) average of all observations π’šn\bm{y}_{n}, each weighted by the responsibility that class kk takes for it.

The E step.

Evidently, the expected sufficient statistics are ⟨XΛ‡kβŸ©π‘ΏΛ‡,𝒀{\mathopen{}\mathclose{{}\left\langle{{\check{X}}_{k}}}\right\rangle_{{\bm{% \check{X}}}{},{\bm{Y}}{}}}, ⟨XΛ‡kβ’π’€βŸ©π‘ΏΛ‡,𝒀{\mathopen{}\mathclose{{}\left\langle{{\check{X}}_{k}{\bm{Y}}}}\right\rangle_{% {\bm{\check{X}}}{},{\bm{Y}}{}}}, and ⟨XΛ‡k⁒𝒀⁒𝒀TβŸ©π‘ΏΛ‡,𝒀{\mathopen{}\mathclose{{}\left\langle{{\check{X}}_{k}{\bm{Y}}{\bm{Y}}^{\text{T% }}}}\right\rangle_{{\bm{\check{X}}}{},{\bm{Y}}{}}}. During EM, i.e.Β when the averaging distribution is p^(𝒙|π’š;𝜽old)p(π’š){\hat{p}\mathopen{}\mathclose{{}\left({\color[rgb]{.75,0,.25}\definecolor[% named]{pgfstrokecolor}{rgb}{.75,0,.25}\bm{x}}{}\middle|{\color[rgb]{.75,0,.25}% \definecolor[named]{pgfstrokecolor}{rgb}{.75,0,.25}\bm{y}};\bm{\theta}^{\text{% old}}}\right)}{p\mathopen{}\mathclose{{}\left({\color[rgb]{.75,0,.25}% \definecolor[named]{pgfstrokecolor}{rgb}{.75,0,.25}\bm{y}}{}}\right)}, these are written more explicitly as

βŸ¨π”ΌXΛ‡k|𝒀[XΛ‡k|𝒀]βŸ©π’€,\displaystyle{\mathopen{}\mathclose{{}\left\langle{\mathbb{E}_{{\check{X}}_{k}% {}|{\bm{Y}}}{\mathopen{}\mathclose{{}\left[{\check{X}}_{k}\middle|{\bm{Y}}{}}% \right]}}}\right\rangle_{{\bm{Y}}{}}},
βŸ¨π”ΌXΛ‡k|𝒀[XΛ‡k|𝒀]π’€βŸ©π’€,\displaystyle{\mathopen{}\mathclose{{}\left\langle{\mathbb{E}_{{\check{X}}_{k}% {}|{\bm{Y}}}{\mathopen{}\mathclose{{}\left[{\check{X}}_{k}\middle|{\bm{Y}}{}}% \right]}{\bm{Y}}}}\right\rangle_{{\bm{Y}}{}}},
βŸ¨π”ΌXΛ‡k|𝒀[XΛ‡k|𝒀]𝒀𝒀TβŸ©π’€.\displaystyle{\mathopen{}\mathclose{{}\left\langle{\mathbb{E}_{{\check{X}}_{k}% {}|{\bm{Y}}}{\mathopen{}\mathclose{{}\left[{\check{X}}_{k}\middle|{\bm{Y}}{}}% \right]}{\bm{Y}}{\bm{Y}}^{\text{T}}}}\right\rangle_{{\bm{Y}}{}}}.

We derived the posterior mean that occurs in all of these expressions in Section 3.1.1. We repeat Eqs.Β 3.7 and 3.8 here for convenience:

𝔼XΛ‡k|𝒀[XΛ‡k|π’š]=softmax{𝒛}k,zk=logΟ€k+12log|𝚺kβˆ’1|βˆ’12(π’šβˆ’πk)T𝚺kβˆ’1(π’šβˆ’πk).\mathbb{E}_{{\check{X}}_{k}{}|{\bm{Y}}}{\mathopen{}\mathclose{{}\left[{\check{% X}}_{k}\middle|\bm{y}{}}\right]}=\operatorname*{softmax}\mathopen{}\mathclose{% {}\left\{\bm{z}}\right\}_{k},\qquad z_{k}=\log\pi_{k}+\frac{1}{2}\log|\mathbf{% {\Sigma}}^{-1}_{k}|-\frac{1}{2}\mathopen{}\mathclose{{}\left(\bm{y}-\bm{\mu}_{% k}}\right)^{\text{T}}\mathbf{{\Sigma}}^{-1}_{k}\mathopen{}\mathclose{{}\left(% \bm{y}-\bm{\mu}_{k}}\right).

9.1.1 KK-means

In Section 3.1.1, we saw what happens to the posterior of the GMM when all classes use the same covariance matrix, 𝚺\mathbf{{\Sigma}}. The class boundaries become lines, and the responsibility of class jj for observation π’š\bm{y} becomes

equation (9.1) (9.1)
p^(X^j=1|π’š;𝜽)=exp⁑{βˆ’12⁒(π’šβˆ’πj)Tβ’πšΊβˆ’1⁒(π’šβˆ’πj)}⁒πjβˆ‘k=1Kexp⁑{βˆ’12⁒(π’šβˆ’πk)Tβ’πšΊβˆ’1⁒(π’šβˆ’πk)}⁒πk.\begin{split}{\hat{p}\mathopen{}\mathclose{{}\left({\hat{X}}_{j}=1\middle|\bm{% y};\bm{\theta}}\right)}&{}=\frac{\exp\mathopen{}\mathclose{{}\left\{-\frac{1}{% 2}\mathopen{}\mathclose{{}\left(\bm{y}-\bm{\mu}_{j}}\right)^{\text{T}}\mathbf{% {\Sigma}}^{-1}\mathopen{}\mathclose{{}\left(\bm{y}-\bm{\mu}_{j}}\right)}\right% \}\pi_{j}}{\sum_{k=1}^{{K}}\exp\mathopen{}\mathclose{{}\left\{-\frac{1}{2}% \mathopen{}\mathclose{{}\left(\bm{y}-\bm{\mu}_{k}}\right)^{\text{T}}\mathbf{{% \Sigma}}^{-1}\mathopen{}\mathclose{{}\left(\bm{y}-\bm{\mu}_{k}}\right)}\right% \}\pi_{k}}.\end{split}

It is not hard to see that in the limit of infinite precision, this quantity goes to zero unless π’š\bm{y} is closer to 𝝁j\bm{\mu}_{j} than any other mean, in which case it goes to 1. The prior probabilities 𝝅\bm{\pi} have become irrelevant. Then the algorithm becomes

KK-MeansΒ Β Β Β Β Β Β Β Β Β Β Β Β Β Β Β Β Β Β Β 

βˆ™\bullet\> E step: pΛ‡(i+1)(X^j=1|π’š)←{1,if ⁒j=argminkβˆ₯π’šβˆ’πk(i)βˆ₯0,otherwise{\check{p}^{(i+1)}\mathopen{}\mathclose{{}\left({\hat{X}}_{j}=1\middle|\bm{y}}% \right)}\leftarrow\begin{cases}1,&\text{if }j=\operatorname*{argmin}_{k}% \mathopen{}\mathclose{{}\left\lVert\bm{y}-\bm{\mu}_{k}^{(i)}}\right\rVert\\ 0,&\text{otherwise}\end{cases}
βˆ™\bullet\> M step: 𝝁j(i)β†βŸ¨pΛ‡(i+1)(X^j=1|𝒀)π’€βŸ©π’€βŸ¨pΛ‡(i+1)(X^j=1|𝒀)βŸ©π’€\bm{\mu}_{j}^{(i)}\leftarrow\frac{{\mathopen{}\mathclose{{}\left\langle{{% \check{p}^{(i+1)}\mathopen{}\mathclose{{}\left({\hat{X}}_{j}=1\middle|{\bm{Y}}% }\right)}{\bm{Y}}}}\right\rangle_{{\bm{Y}}{}}}}{{\mathopen{}\mathclose{{}\left% \langle{{\check{p}^{(i+1)}\mathopen{}\mathclose{{}\left({\hat{X}}_{j}=1\middle% |{\bm{Y}}}\right)}}}\right\rangle_{{\bm{Y}}{}}}}

This algorithm, which pre-existed EM for the GMM, is known as KK-means††margin: KK-means .