高斯混合模型(Gaussian Mixture Model)與EM演算法原理(二)

接前文。

6. 高斯混合模型的EM演算法

前文講到，不完全數據集的對數似然函數 $ln p(mathbf{X|pi,mu,Sigma}) = sum_{n=1}^{N}{lnleft[ sum_{k=1}^{K}{pi_k pleft(mathbf{x_n | mu_k, Sigma_k} ight)} ight]}$ 的極大似然估計得到參數：

$left{ egin{align} mathbf{mu_k} &= frac{1}{N_k} sum_{n=1}^{N}{left[ gamma(z_{nk})mathbf{x_n} ight]}\ mathbf{Sigma_k} &= frac{1}{N_k} sum_{n=1}^{N}{gamma(z_{nk})(mathbf{x_n-mu_k})(mathbf{x_n-mu_k})^T} \ pi_k &= frac{N_k}{N} end{align}\ ight. ag{1}$

其中，響應度 $gamma(z_{nk}) = frac{pi_k pleft(mathbf{x_n | mu_k, Sigma_k} ight)} {sum_{j=1}^{K}{pi_j pleft(mathbf{x_n | mu_j, Sigma_j} ight)} }$ ， $N_k = sum_{n=1}^{N}{gamma(z_{nk})}$ 。

雖然(1)式不是閉式解，但是我們可以根據(1)式，通過迭代的方法來計算參數 $mathbf{mu_k, Sigma_k}, pi_k$ 。

EM演算法(Expectation-Maximization algorithm)是就是這樣一種用於計算含有隱變數的模型的極大似然解的強大方法。

下面簡述用EM演算法來計算高斯混合模型參數的步驟：

1. 初始化：給參數均值向量 $mathbf{mu_k}$ ，協方差矩陣 $mathbf{Sigma_k}$ 和混合係數賦初值；並計算對數似然函數的初值；

* 註：可以用K均值聚類(K-means clustering)演算法來得到初始參數。

2. E步(Expectation step)：用當前參數值 $mathbf{mu_k, Sigma_k}, pi_k$ 計算後驗概率(響應度) $gamma(z_{nk})$ ：

$gamma(z_{nk}) leftarrow frac{pi_k pleft(mathbf{x_n | mu_k, Sigma_k} ight)} {sum_{j=1}^{K}{pi_j pleft(mathbf{x_n | mu_j, Sigma_j} ight)} }$

3. M步(Maximization step)：用當前響應度 $gamma(z_{nk})$ ，根據(1)式重新估計參數：

$egin{align} mathbf{mu_k^{new}} &leftarrow frac{1}{N_k} sum_{n=1}^{N}{left[ gamma(z_{nk})mathbf{x_n} ight]}\ mathbf{Sigma_k^{new}} &leftarrow frac{1}{N_k} sum_{n=1}^{N}{gamma(z_{nk})(mathbf{x_n-mu_k^{new}})(mathbf{x_n-mu_k^{new}})^T} \ pi_k^{new} &leftarrow frac{N_k}{N} end{align}$

其中， $N_k = sum_{n=1}^{N}{gamma(z_{nk})}$ 。

在M步中，先計算 $mathbf{mu_k^{new}}$ 的值，然後用 $mathbf{mu_k^{new}}$ 來計算新的協方差矩陣 $mathbf{Sigma_k^{new}}$ 。

4. 用新的參數 $mathbf{mu_k^{new}}, mathbf{Sigma_k^{new}}, pi_k^{new}$ 計算對數似然函數：

$ln p(mathbf{X|pi,mu,Sigma}) = sum_{n=1}^{N}{lnleft[ sum_{k=1}^{K}{pi_k pleft(mathbf{x_n | mu_k, Sigma_k} ight)} ight]}$

然後檢查是否收斂。(當然，也可直接檢查參數是否收斂。)

如果收斂，則得到了模型參數；如果還沒收斂，則返回第2步繼續迭代。

參考文獻[1], Section9.2, Page437指出，E步+M步可以保證對數似然函數的上升。不過需要注意的是，作為一種迭代演算法，EM演算法有可能收斂到局部最大值(local maximum)。

7. 換個角度看高斯混合模型的EM演算法

在前面的極大似然估計中，我們的目標函數是不完全數據集 $mathbf{X=left{ x_1, ldots, x_N ight}}$ 的對數似然函數。

現在我們換個角度來看：假如我們也有隱變數的觀測值 $mathbf{Z=left{ z_1, ldots, z_N ight}}$ ，即知道 $mathbf{x_n}$ 來自高斯混合中的哪一個分模型，換句話說，我們有完全數據集，那麼我們可以計算完全數據的對數似然函數。

首先，一組完全數據 $left{ mathbf{x_n, z_n} ight}$ 的似然函數為：

$p(mathbf{x_n, z_n|pi,mu,Sigma}) = prod_{k=1}^{K}{pi_k^{z_{nk}} pleft(mathbf{x_n | mu_k, Sigma_k} ight)^{z_{nk}}}$

所以，整個完全數據集的似然函數為：

$p(mathbf{X, Z|pi,mu,Sigma}) = prod_{n=1}^{N}prod_{k=1}^{K}{pi_k^{z_{nk}} pleft(mathbf{x_n | mu_k, Sigma_k} ight)^{z_{nk}}}$

(當然，我們假設每組數據相互獨立。)

取對數得

$ln p(mathbf{X, Z|pi,mu,Sigma}) = sum_{n=1}^{N}{ sum_{k=1}^{K}{left{ {z_{nk}left[ln pi_k + ln pleft(mathbf{x_n | mu_k, Sigma_k} ight) ight]} ight}} } ag{2}$

可以看出，和不完全數據的對數似然函數 $ln p(mathbf{X|pi,mu,Sigma}) = sum_{n=1}^{N}{lnleft[ sum_{k=1}^{K}{pi_k pleft(mathbf{x_n | mu_k, Sigma_k} ight)} ight]}$ 相比，(2)式中 $sum_{k=1}^{K}$ 和換了順序。

然而實際情況是，我們並不知道隱變數的觀測值，即(2)式中的 $z_{nk}$ ，所以我們無法直接計算。因此，我們考慮對在隱變數 $z_{nk}$ 的後驗分布下求期望，利用求期望的過程，將 $z_{nk}$ 消掉。

在隱變數 $z_{nk}$ 的後驗分布下的期望：

$egin{align} E_{mathbf Z} ln p(mathbf{X, Z|pi,mu,Sigma}) &= E sum_{n=1}^{N}{ sum_{k=1}^{K}{left{ {z_{nk}left[ln pi_k + ln pleft(mathbf{x_n | mu_k, Sigma_k} ight) ight]} ight}}}\ &=sum_{n=1}^{N}{ sum_{k=1}^{K}{left{E(z_{nk})left[ln pi_k + ln pleft(mathbf{x_n | mu_k, Sigma_k} ight) ight] ight}}} end{align}$

其中， $E(z_{nk})$ 表示 $z_{nk}$ 在後驗分布下的期望 $E(z_{nk}|mathbf{X, heta})$ 。因 $z_{nk}in {0, 1}$ ，所以

$egin{align} E(z_{nk}|mathbf{X, heta}) &= 0 imes p(z_{nk}=0|mathbf{X, heta}) + 1 imes p(z_{nk}=1|mathbf{X, heta})\ &= p(z_{nk}=1|mathbf{X, heta})\ &= gamma(z_{nk}) end{align}$

也就是說，期望 $E(z_{nk})$ 就是響應度 $gamma(z_{nk})$ 。進而，

$egin{align} E_{mathbf Z}ln p(mathbf{X, Z|pi,mu,Sigma}) &=sum_{n=1}^{N}{ sum_{k=1}^{K}{left{gamma(z_{nk})left[ln pi_k + ln pleft(mathbf{x_n | mu_k, Sigma_k} ight) ight] ight}}} end{align}$

將上式簡記為 $E_{mathbf{Z}}$ ，現在將 $gamma(z_{nk})$ 當作常數，求參數，使 $E_{mathbf{Z}}$ 最大化，即

$max_{mathbf{pi,mu,Sigma}}{E_{mathbf{Z}}}\ s.t. sum_{k=1}^{K}{pi_k} = 1$

採用拉格朗日乘子法，拉格朗日函數為：

$L_E = E + lambda left( sum_{k=1}^{K}{pi_k} - 1 ight)$

先對 $mathbf{mu_k}$ 求梯度，因為

$egin{align} frac{partial ln pleft(mathbf{x_n | mu_k, Sigma_k} ight)}{partial mathbf{mu_k}} &= frac{1}{pleft(mathbf{x_n | mu_k, Sigma_k} ight)} pleft(mathbf{x_n | mu_k, Sigma_k} ight)left[ mathbf{Sigma^{-1}(x_n-mu_k)} ight] \&= mathbf{Sigma^{-1}(x_n-mu_k)} end{align}$

所以，

$egin{align} frac{partial L_E}{partial mathbf{mu_k}} &= sum_{n=1}^{N}{ left[ frac{partial sum_{k=1}^{K}{gamma(z_{nk})ln pleft(mathbf{x_n | mu_k, Sigma_k} ight)}}{partial mathbf{mu_k}} ight]}\ &= sum_{n=1}^{N}{ left[gamma(z_{nk}) frac{partial ln pleft(mathbf{x_n | mu_k, Sigma_k} ight)}{partial mathbf{mu_k}} ight]}\ &= sum_{n=1}^{N}{gamma(z_{nk})left[ mathbf{Sigma_k^{-1}(x_n-mu_k)} ight]}\ end{align}$

進而，

$egin{align} frac{partial L_E}{partial mathbf{mu_k}} = mathbf{0} &Rightarrow sum_{n=1}^{N}{gamma(z_{nk})left[ mathbf{Sigma_k^{-1}(x_n-mu_k)} ight]} = mathbf{0}\ &Rightarrow sum_{n=1}^{N}{gamma(z_{nk})mathbf{(x_n-mu_k)}} = mathbf{0}\ &Rightarrow sum_{n=1}^{N}{gamma(z_{nk})mathbf{x_n}} = sum_{n=1}^{N}{gamma(z_{nk})mathbf{mu_k}}\ end{align}$

同之前極大化一樣，可推得

$mathbf{mu_k} = frac{1}{N_k} sum_{n=1}^{N}{left[ gamma(z_{nk})mathbf{x_n} ight]}$

然後是對 $mathbf{Sigma_k}$ 求梯度，

$egin{align} frac{partial L_E}{partial mathbf{Sigma_k}} &= sum_{n=1}^{N}{ left[ frac{partial sum_{k=1}^{K}{gamma(z_{nk})ln pleft(mathbf{x_n | mu_k, Sigma_k} ight)}}{partial mathbf{Sigma_k}} ight]}\ &= sum_{n=1}^{N}{ left[gamma(z_{nk}) frac{partial ln pleft(mathbf{x_n | mu_k, Sigma_k} ight)}{partial mathbf{Sigma_k}} ight]}\ &= sum_{n=1}^{N}{ left[frac{gamma(z_{nk})}{pleft(mathbf{x_n | mu_k, Sigma_k} ight)} frac{partial pleft(mathbf{x_n | mu_k, Sigma_k} ight)}{partial mathbf{Sigma_k}} ight]}\ end{align} ag{3}$

回憶在文(一)中，我們求的極值點時，有

$egin{align} frac{partial L}{partial mathbf{Sigma_k}} = sum_{n=1}^{N}{left[ frac{pi_k }{sum_{j=1}^{K}{pi_j pleft(mathbf{x_n | mu_j, Sigma_j} ight)}} frac{partial pleft(mathbf{x_n | mu_k, Sigma_k} ight)}{partial mathbf{Sigma_k}} ight]} end{align} ag{4}$

對比(3)(4)兩式，結合響應度 $gamma(z_{nk})$ 的公式，可發現兩式是一致的，這意味著兩種方法計算出的 $mathbf{Sigma_k}$ 也相同。

再來看，

$egin{align} frac{partial L_E}{partial pi_k} &= sum_{n=1}^{N}{ left[ frac{partial sum_{k=1}^{K}{gamma(z_{nk})ln pi_k}}{partial pi_k} ight]} + lambda\ &= sum_{n=1}^{N}{ left[gamma(z_{nk}) frac{partial ln pi_k}{partial pi_k} ight]} + lambda\ &= sum_{n=1}^{N}{ left[gamma(z_{nk}) frac{1}{pi_k} ight]} + lambda = frac{1}{pi_k}sum_{n=1}^{N}{gamma(z_{nk})} + lambda\ &= frac{1}{pi_k} N_k + lambda end{align}$

對比前文，可以發現也和極值點的參數一致。

綜上，通過極大化 $E_{mathbf Z}ln p(mathbf{X, Z|pi,mu,Sigma})$ 得到的模型參數，和的極大似然估計結果是一致的。在EM演算法的每一次迭代中，我們都在對 $E_{mathbf{Z}}$ 求極大值。

這給了我們另一個角度來看EM演算法在高斯混合模型中的應用，也有助於我們理解EM演算法的一般形式。

8. EM演算法的一般形式

對於含有隱變數的概率模型，記變數和(離散的)隱變數分別為，模型參數為，(比如高斯混合模型的參數 )， $mathbf{X = {x_1,ldots, x_N}}$ 為變數的個觀測值的集合，即不完全數據集， $mathbf{Z = {z_1,ldots, z_N}}$ 為相應的隱變數值的集合。

要計算極大化對數似然函數，EM演算法的思路是：考慮完全數據集的對數似然函數在隱變數的後驗分布 $p(mathbf{Z | X, heta^{old}})$ 下的期望，反覆迭代極大化這個期望，最終得到模型的參數值。其中 $mathbf{ heta^{old}}$ 為當前參數。(若 $mathbf{ x_1,ldots, x_N}$ 相互獨立，則計算 $p(mathbf{Z | X, heta^{old}})$ 實際上是計算 $p(mathbf{z_n | x_n, heta^{old}}), n = 1,ldots, N$ )