前言
在模擬不明隨機變數時通常會選用 Gaussian 作為該變數的分布，因為科學家發現各式各樣的物理現象都神奇地符合 Gaussian 分布，但是生活中也會遇到非單一鐘形分布的變數，例如圖像的色彩分布、語音的發聲分布等，這些分布比較像多個鐘形分布的結合，因此 Gaussian Mixture Model (GMM) 被提出，為了可以擬合單一 Gaussian distribution 無法擬合的隨機變數，本文章將介紹 GMM 模型，以及如何利用 EM algorithm 求出 GMM 的參數，內含詳細的推導過程，需有基礎的線性代數運算能力才能夠融會貫通。
Gaussian Mixture Model
GMM 模型是多個 Gaussian distribution 的線性組合，其數學式如下：$$
p(\textbf{x}) = \sum^{K}_{k=1} \pi_k \mathcal{N}(\textbf{x}|\mathbf{\mu}_k, \mathbf{\Sigma}_k)
$$ 代表此 GMM 模型由 $K$ 個 Gaussian distribution 結合而成。為了使 GMM 符合機率分布要求，即積分為 $1$，Gaussian 的權重 $\pi_k$ 總和必須為 $1$: $$
\sum^{K}_{k=1} \pi_k = 1
$$ 讀者可以利用簡易的積分檢查此 GMM 是否符合積分為 $1$ 之要求。這裡我們引入另一個 one-hot encoded 的隨機向量 $\textbf{z}$，且該向量的機率分布為：$$
p(\textbf{z}) = \prod^K_{k=1} \pi_k \Rightarrow p(z_k = 1) = \pi_k
$$ 因此我們可以重寫 GMM 模型：$$
p(\textbf{x}) = \sum^{K}_{k=1} \pi_k \mathcal{N}(\textbf{x}|\mathbf{\mu}_k, \mathbf{\Sigma}_k) = \sum_{\textbf{z}} p(\textbf{z}) p(\textbf{x}|\textbf{z}) = \sum_{\textbf{z}} p(\textbf{x}, \textbf{z})
$$ 此 $\pi_k$ 可以被看作 $\textbf{x}$ 來自於第 $k$ 群的機率，且 $\mathcal{N}(\textbf{x}|\mathbf{\mu}_k, \mathbf{\Sigma}_k)$ 代表屬於第 $k$ 群的條件下 $\textbf{x}$ 的機率分布。這裡我雖然用”屬於”這個詞，但是實際上 $\textbf{x}$ 可能包含各群的成分在裡面，不一定非黑即白地屬於某個群。
EM Algorithm
要求得 GMM 的變數 $\mathbf{\mu}_k$、$\mathbf{\Sigma}_k$ 困難之處在於，我們無法直接利用 ML estimation 算出封閉解，理由是我們只能夠蒐集到 $\{\textbf{x}_1, …, \textbf{x}_N\}$ 的資料，卻無法蒐集到 $\{\textbf{z}_1, …, \textbf{z}_N\}$，因為 $\textbf{z}_k$ 是一個 latent variable，我們無從得知該變數的實際值為何。為了求出 GMM 的變數，我們必須從 $p(\textbf{x}, \textbf{z})$ 出發並利用 EM algorithm 估計，就如同丟硬幣的例子，利用 posterior probability 估計出 maximum likelihood 的值：$$
\mathbb{E} \bigg[ \ln p(\textbf{X}, \textbf{Z})\bigg] = \sum_{i=1}^{N} \sum_{\textbf{z}_i} p(\textbf{z}_i | \textbf{x}_i) p(\textbf{x}_i, \textbf{z}_i)
$$ 其中我們可以利用以下式子得出 $p(\textbf{z}_i | \textbf{x}_i)$ 的數學式：$$
\begin{aligned}
p(\textbf{z}_i | \textbf{x}_i) &= \frac{p(\textbf{z}_i)p(\textbf{x}_i | \textbf{z}_i)}{\sum_{\textbf{z}_i} p(\textbf{z}_i)p(\textbf{x}_i | \textbf{z}_i)} \\
&= \frac{\pi_k \mathcal{N}(\textbf{x}_i | \mathbf{\mu}_k, \mathbf{\sigma}_k)}{\sum_{j=1}^K \pi_j \mathcal{N}(\textbf{x}_i | \mathbf{\mu}_j, \mathbf{\sigma}_j)}
\end{aligned}$$ 因此 ML 的估計值可以重寫為：$$
\mathbb{E} \bigg[ \ln p(\textbf{X}, \textbf{Z})\bigg] = \sum_{i=1}^{N} \sum_{k=1}^{K} \frac{\pi_k \mathcal{N}(\textbf{x}_i | \mathbf{\mu}_k, \mathbf{\sigma}_k)}{\sum_{j=1}^K \pi_j \mathcal{N}(\textbf{x}_i | \mathbf{\mu}_j, \mathbf{\sigma}_j)} \ln \pi_k \mathcal{N}(\textbf{x}_i|\mathbf{\mu}_k, \mathbf{\Sigma}_k) $$
E step
E step 主要目的是算出 posterior probability 的實際值並且固定，以進行下一步的估計，為了得出此機率值，我們必須有參數的起始值，因此我們會先假設 $\pi_k, \mathbf{\mu}_k, \mathbf{\Sigma}_k, k \in {1, …, K}$ 的初始值並且算出 posterior probability，我們以 $\gamma_{ik}$ 表示：$$
\gamma_{ik} = \frac{\pi_k \mathcal{N}(\textbf{x}_i | \mathbf{\mu}_k, \mathbf{\sigma}_k)}{\sum_{j=1}^K \pi_j \mathcal{N}(\textbf{x}_i | \mathbf{\mu}_j, \mathbf{\sigma}_j)}
$$
M step – $\pi_k$
計算出 posterior probability $\gamma_{ik}$ 值後，我們可以開始利用偏微分分別算出 $\pi_k, \mathbf{\mu}_k, \mathbf{\Sigma}_k$ 的估計值。首先我們對 $\pi_k$ 做偏微分並設為 $0$，這裡必須注意的是此最佳化是一個受限制的最佳化，因為 $\sum_{k=1}^K \pi_k = 1$，因此我們需要利用 Lagrange multiplier 幫助我們求解：$$
\begin{aligned}
&\frac{\partial}{\partial \pi_k} \bigg\{ \sum_{i=1}^{N} \sum_{t=1}^{K} \gamma_{it} \ln \pi_t \mathcal{N}(\textbf{x}_i|\mathbf{\mu}_t, \mathbf{\Sigma}_t) + \lambda \sum_{t=1}^K \pi_t \bigg\} = 0\\
\Rightarrow \hspace{1cm} &\frac{1}{\pi_k} \sum_{i=1}^{N} \gamma_{ik} + \lambda = 0 \\
\Rightarrow \hspace{1cm} &\pi_k = \frac{\sum_{i=1}^{N} \gamma_{ik}}{-\lambda} \\
\Rightarrow \hspace{1cm} &\sum_{k=1}^{K} \pi_k = \sum_{k=1}^{K} \frac{\sum_{i=1}^{N} \gamma_{ik}}{-\lambda} = \frac{N}{-\lambda} = 1 \Rightarrow N = -\lambda \\
\Rightarrow \hspace{1cm} &\pi_k = \frac{\sum_{i=1}^{N} \gamma_{ik}}{N}
\end{aligned}
$$ 利用 posterior probability 的總和為 1: $\sum_{k=1}^K \gamma_{ik} = 1$，即可求得第四個等式的結果。
M step – $\mathbf{\mu}_k$
我們對 $\mathbf{\mu}_k$ 做偏微分之前，我們必須知道以下等式：$$
\frac{\partial \mathbf{y}^\top \mathbf{B} \mathbf{y}}{\partial \mathbf{y}} = (\mathbf{B} + \mathbf{B}^\top) \mathbf{y}, \hspace{0.5cm} \text{if} \hspace{0.5cm} \mathbf{B}^\top = \mathbf{B} \hspace{0.5cm} \text{then} \hspace{0.5cm} \frac{\partial \mathbf{y}^\top \mathbf{B} \mathbf{y}}{\partial \mathbf{y}} = 2\mathbf{B} \mathbf{y}
$$ 藉由此向量微分等式，我們可以算出對 $\mathbf{\mu}_k$ 之偏微分：$$
\begin{aligned}
&\frac{\partial}{\partial \mathbf{\mu}_k} \sum_{i=1}^{N} \sum_{t=1}^{K} \gamma_{it} \ln \pi_t \mathcal{N}(\textbf{x}_i|\mathbf{\mu}_t, \mathbf{\Sigma}_t) = 0 \\
\Rightarrow \hspace{1cm} &\frac{\partial}{\partial \mathbf{\mu}_k} \sum_{i=1}^{N} \sum_{t=1}^{K} \gamma_{tk} \ln \frac{\pi_t}{\sqrt{(2 \pi)^n |\mathbf{\Sigma}_t|}} \exp \bigg\{ \frac{-1}{2} (\mathbf{x}_i – \mathbf{\mu}_t)^\top \mathbf{\Sigma}_t^{-1} (\mathbf{x}_i – \mathbf{\mu}_t) \bigg\}\\
\Rightarrow \hspace{1cm} &\sum_{i=1}^{N} \gamma_{ik} \frac{\pi_k \mathcal{N}(\textbf{x}|\mathbf{\mu}_k, \mathbf{\Sigma}_k)}{\pi_k \mathcal{N}(\textbf{x}|\mathbf{\mu}_k, \mathbf{\Sigma}_k)} \cdot (-\frac{1}{2}) \cdot 2 \cdot \mathbf{\Sigma}_k^{-1} (\mathbf{x}_i – \mathbf{\mu}_k) (-1) \\
\Rightarrow \hspace{1cm} &\sum_{i=1}^{N} \gamma_{ik} \mathbf{\Sigma}_k^{-1} (\mathbf{x}_i – \mathbf{\mu}_k) = 0 \\
\Rightarrow \hspace{1cm} &\mathbf{\mu}_k = \frac{\sum_{i=1}^{N} \gamma_{ik} \mathbf{x}_i}{\sum_{i=1}^{N} \gamma_{ik}}
\end{aligned}
$$ 其中 covariance matrix $\mathbf{\Sigma}$ 是一個 symmetric matrix，所以 $\mathbf{\Sigma}^{-1}$ 也是一個 symmetric matrix，第一個等式到第二個等式推倒過程中利用了這個性質，並搭配簡單的微分計算即可求得。
M step – $\mathbf{\Sigma}_k$
最後我們對 $\mathbf{\Sigma}$ 做偏微分可以求出其估計值，在此之前我們先介紹幾個等式：$$
\begin{aligned}
&1: \frac{\partial \mathbf{a}^\top \mathbf{X}^{-1} \mathbf{b}}{\partial \mathbf{X}} = -\mathbf{X}^{-\top} \mathbf{a} \mathbf{b}^\top \mathbf{X}^{-\top} \\
&2: \mathbf{X}^{-1} = \frac{\text{adj} (\mathbf{X})}{|\mathbf{X}|}, \hspace{0.5cm} \text{if} \hspace{0.5cm} \mathbf{X} = \mathbf{X}^\top \hspace{0.5cm} \text{then} \hspace{0.5cm} \mathbf{X}^{-1} = \frac{\text{adj}^\top (\mathbf{X})}{|\mathbf{X}|} \\
&3: \frac{d |\mathbf{X}|}{d \mathbf{X}} = \text{adj}^\top (\mathbf{X})
\end{aligned}
$$ 如同前面所說， $\mathbf{\Sigma}$ 是對稱矩陣，因此等式 $2$ 中的 $\mathbf{\Sigma}^\top = \mathbf{\Sigma}$ 成立。至此我們可以開始推導 $\mathbf{\Sigma}$ 的偏微分了：$$
\begin{aligned}
&\frac{\partial}{\partial \mathbf{\Sigma}_k} \sum_{i=1}^{N} \sum_{t=1}^{K} \gamma_{tk} \ln \pi_t \mathcal{N}(\textbf{x}_i|\mathbf{\mu}_t, \mathbf{\Sigma}_t) = 0 \\
\Rightarrow \hspace{1cm} &\frac{\partial}{\partial \mathbf{\Sigma}_k} \sum_{i=1}^{N} \sum_{t=1}^{K} \gamma_{ik} \ln \frac{\pi_t}{\sqrt{(2 \pi)^n |\mathbf{\Sigma}_t|}} \exp \bigg\{ \frac{-1}{2} (\mathbf{x}_i – \mathbf{\mu}_t)^\top \mathbf{\Sigma}_t^{-1} (\mathbf{x}_i – \mathbf{\mu}_t) \bigg\} = 0\\
\Rightarrow \hspace{1cm} &\frac{\partial}{\partial \mathbf{\Sigma}_k} \bigg\{ \sum_{i=0}^N \gamma_{ik} \ln \frac{\pi_k}{\sqrt{(2 \pi)^n}} -\frac{1}{2} \sum_{i=0}^N \gamma_{ik} \ln |\mathbf{\Sigma}_k| -\frac{1}{2} \sum_{i=0}^N \gamma_{ik} (\mathbf{x}_i – \mathbf{\mu}_k)^\top \mathbf{\Sigma}_k^{-1} (\mathbf{x}_i – \mathbf{\mu}_k) \bigg\} = 0 \\
\Rightarrow \hspace{1cm} &-\frac{1}{2} \sum_{i=0}^N \gamma_{ik} \frac{\text{adj} (\mathbf{\Sigma}_k)}{|\mathbf{\Sigma}_k|} + \frac{1}{2} \sum_{i=0}^N \gamma_{ik} \mathbf{\Sigma}_k^{-\top} (\mathbf{x}_i – \mathbf{\mu}_k) (\mathbf{x}_i – \mathbf{\mu}_k)^\top \mathbf{\Sigma}_k^{-\top} = 0 \\
\Rightarrow \hspace{1cm} &-\frac{1}{2} \sum_{i=0}^N \gamma_{ik} \mathbf{\Sigma}_k^{-1} + \frac{1}{2} \sum_{i=0}^N \gamma_{ik} \mathbf{\Sigma}_k^{-T} (\mathbf{x}_i – \mathbf{\mu}_k) (\mathbf{x}_i – \mathbf{\mu}_k)^\top \mathbf{\Sigma}_k^{-T} = 0 \\
\Rightarrow \hspace{1cm} &\sum_{i=0}^N \gamma_{ik} \mathbf{\Sigma}_k^{-1} (\mathbf{x}_i – \mathbf{\mu}_k) (\mathbf{x}_i – \mathbf{\mu}_k)^\top \mathbf{\Sigma}_k^{-1} = \sum_{i=0}^N \gamma_{ik} \mathbf{\Sigma}_k^{-1} \\
\Rightarrow \hspace{1cm} &\mathbf{\Sigma}_k = \frac{\sum_{i=0}^N \gamma_{ik} (\mathbf{x}_i – \mathbf{\mu}_k) (\mathbf{x}_i – \mathbf{\mu}_k)^\top}{\sum_{i=0}^N \gamma_{ik}}
\end{aligned}
$$ 至此我們已經把所有參數的迭代算式寫出，剩下的只要依照 EM algorithm 的步驟慢慢逼近局部最佳解即可。
其他文章
機器學習系列文章總整理
Maximum Likelihood & Maximum a Posteriori－基礎估計模型的詳細介紹
EM Algorithm 詳盡介紹: 利用簡單例子輕鬆讀懂 EM 的原理及概念
[MobileNet 系列] MobileNetV1－為終端而生
[MobileNet 系列] MobileNetV2－小改動大進化，終端的唯一選擇
Pingback: 機器學習系列文章總整理 - PlayRound 玩轉部落格