EM算法及其應用（一）

時間 2019-12-14

標籤算法及其應用欄目 CSS 简体版

原文原文鏈接

EM算法及其應用（一）

EM算法及其應用（二）： K-means 與高斯混合模型

EM算法是指望最大化 (Expectation Maximization) 算法的簡稱，用於含有隱變量的狀況下，機率模型參數的極大似然估計或極大後驗估計。EM算法是一種迭代算法，每次迭代由兩步組成：E步，求指望 (expectation)，即利用當前估計的參數值來計算對數似然函數的指望值；M步，求極大 (maximization)，即求參數$\theta$ 來極大化E步中的指望值，而求出的參數$\theta$將繼續用於下一個E步中指望值的估計。EM算法在機器學習中應用普遍，本篇和下篇文章分別探討EM算法的原理和其兩大應用 —— K-means和高斯混合模型。html

$\large{\S} \normalsize\mathrm{1}$ 先驗知識

凸函數、凹函數和 Jensen不等式算法

設$f(x)$爲定義在區間$I = [a,b]$上的實值函數，對於任意$\forall \, x_1, x_2 \in I, \lambda \in [0,1]$，有：
\[ f(\lambda \,x_1 + (1-\lambda)\,x_2) \leq \lambda f(x_1) + (1-\lambda)f(x_2) \]
則$f(x)$爲凸函數 (convex function)，以下圖所示。相應的，若上式中 $\leqslant$ 變爲 $\geqslant$ ，則$f(x)$爲凹函數 (concave function)。凸函數的斷定條件是二階導 $f^{''}(x) \geqslant 0$，而凹函數爲 $f^{''}(x) \leqslant 0$ 。後文要用到的對數函數$ln(x)$的二階導爲$-\frac{1}{x^2} < 0$，因此是凹函數。機器學習

Jensen不等式就是上式的推廣，設$f(x)$爲凸函數，$\lambda_i \geqslant 0, \;\; \sum_i \lambda_i = 1$，則：
\[ f\left(\sum\limits_{i=1}^n \lambda_i x_i\right) \leq \sum\limits_{i=1}^n \lambda_i f(x_i) \]
若是是凹函數，則將不等號反向，若用對數函數來表示，就是：
\[ ln\left(\sum\limits_{i=1}^n \lambda_i x_i\right) \geq \sum\limits_{i=1}^n \lambda_i ln(x_i) \]
若將$\lambda_i$視爲一個機率分佈，則可表示爲指望值的形式，在後文中一樣會引入機率分佈：
\[ f(\mathbb{E}[\mathrm{x}]) \leq \mathbb{E}[f(\mathrm{x})] \]函數

KL散度post

KL散度(Kullback-Leibler divergence) 又稱相對熵 (relative entropy)，主要用於衡量兩個機率分佈p和q的差別，也可理解爲兩個分佈對數差的指望。
\[ \mathbb{KL}(p||q) = \sum_i p(x_i)log \frac{p(x_i)}{q(x_i)}= \mathbb{E}_{\mathrm{x}\sim p}\left[log \frac{p(x)}{q(x)}\right] = \mathbb{E}_{\mathrm{x}\sim p}\left[log\,p(x) - log\,q(x) \right ] \]
KL散度總知足$\mathbb{KL}(p||q) \geqslant 0$，而當且僅當$q=p$時，$\mathbb{KL}(p||q) = 0$ 。通常來講分佈$p(x)$比較複雜，於是但願用比較簡單的$q(x)$去近似$p(x)$，而近似的標準就是KL散度越小越好。學習

KL散度不具有對稱性，即$\mathbb{KL}(p||q) \neq \mathbb{KL}(q||p)$，所以不能做爲一個距離指標。優化

極大似然估計和極大後驗估計spa

極大似然估計 (Maximum likelihood estimation) 是參數估計的經常使用方法，基本思想是在給定樣本集的狀況下，求使得該樣本集出現的「可能性」最大的參數$\theta$。將參數$\theta$視爲未知量，則參數$\theta$對於樣本集X的對數似然函數爲：
\[ L(\theta) = ln \,P(X|\theta) \]
這個函數反映了在觀測結果X已知的條件下，$\theta$的各類值的「似然程度」。這裏是把觀測值X當作結果，把參數$\theta$當作是致使這個結果的緣由。參數$\theta$雖然未知可是有着固定值 (固然這是頻率學派的觀點)，並不是事件或隨機變量，無幾率可言，於是改用「似然(likelihood)" 這個詞。orm

因而經過求導求解使得對數似然函數最大的參數$\theta$，$\theta = \mathop{\arg\max}\limits_{\theta}L(\theta)$，即爲極大似然法。htm

極大後驗估計 (Maximum a posteriori estimation) 是貝葉斯學派的參數估計方法，相比於頻率學派，貝葉斯學派將參數$\theta$視爲隨機變量，並將其先驗分佈$P(\theta)$包含在估計過程當中。運用貝葉斯定理，參數$\theta$的後驗分佈爲：
\[ P(\theta|X) = \frac{P(X,\theta)}{P(X)} = \frac{P(\theta)P(X|\theta)}{P(X)} \propto P(\theta)P(X|\theta) \]
上式中$P(X)$不依賴於$\theta$於是爲常數項能夠捨去，則最終結果爲 $\theta = \mathop{\arg\max}\limits_{\theta}P(\theta)P(X|\theta)$

$\large{\S} \normalsize\mathrm{2}$ EM算法初探

機率模型有時既含有觀測變量 (observable variable)，又含有隱變量 (hidden variable)，隱變量顧名思義就是沒法被觀測到的變量。若是都是觀測變量，則給定數據，能夠直接使用極大似然估計。但若是模型含有隱變量時，直接求導獲得參數比較困難。而EM算法就是解決此類問題的經常使用方法。

對於一個含有隱變量$\mathbf{Z}$的機率模型，通常將$\{\mathbf{X}, \mathbf{Z}\}$稱爲徹底數據，而觀測數據$\mathbf{X}$爲不徹底數據。

咱們的目標是極大化觀測數據$\mathbf{X}$關於參數$\boldsymbol{\theta}$的對數似然函數。因爲存在隱變量，於是也可表示爲極大化$\mathbf{X}$的邊緣分佈 (marginal distribution)，即：
\[ L(\boldsymbol{\theta}) = ln\,P(\mathbf{X}|\boldsymbol{\theta}) = ln\,\sum\limits_{\mathbf{Z}}P(\mathbf{X},\mathbf{Z}|\boldsymbol{\theta}) \tag{1.1} \]
上式中存在「對數的和」 —— $ln\sum(\cdot)$，若是直接求導將會很是困難。於是EM算法採用曲線救國的策略，構建$(1.1)$式的一個下界，而後經過極大化這個下界來間接達到極大化$(1.1)$的效果。

要想構建下界，就須要運用上文中的Jensen不等式。記$\boldsymbol{\theta}^{(t)}$爲第t步迭代參數的估計值，考慮引入一個分佈$P(\mathbf{Z}|\mathbf{X},\boldsymbol{\theta}^{(t)})$，因爲：

$P(\mathbf{Z}|\mathbf{X},\boldsymbol{\theta}^{(t)}) \geqslant 0$
$\sum_{\mathbf{Z}}P(\mathbf{Z}|\mathbf{X},\boldsymbol{\theta}^{(t)}) = 1$
$ln(\cdot)$爲凹函數

於是能夠利用Jensen不等式求出$L(\boldsymbol{\theta})$的下界：
\[ \begin{align} L(\boldsymbol{\theta}) = ln\,\sum\limits_{\mathbf{Z}}P(\mathbf{X},\mathbf{Z}|\boldsymbol{\theta}) &= ln\,\sum\limits_{\mathbf{Z}}P(\mathbf{Z|\mathbf{X},\boldsymbol{\theta}^{(t)}})\frac{P(\mathbf{X},\mathbf{Z}|\boldsymbol{\theta}) }{P(\mathbf{Z}|\mathbf{X},\boldsymbol{\theta}^{(t)})} \tag{1.2}\\ & \geqslant \sum\limits_{\mathbf{Z}}P(\mathbf{Z|\mathbf{X},\boldsymbol{\theta}^{(t)}}) \,ln\frac{P(\mathbf{X},\mathbf{Z}|\boldsymbol{\theta}) }{P(\mathbf{Z}|\mathbf{X},\boldsymbol{\theta}^{(t)})} \tag{1.3} \\ & = \underbrace{\sum\limits_{\mathbf{Z}}P(\mathbf{Z|\mathbf{X},\boldsymbol{\theta}^{(t)}}) \,ln{P(\mathbf{X},\mathbf{Z}|\boldsymbol{\theta}) }}_{\mathcal{Q}(\boldsymbol{\theta},\boldsymbol{\theta}^{(t)})} \;\;\underbrace{- \sum\limits_{\mathbf{Z}}P(\mathbf{Z|\mathbf{X},\boldsymbol{\theta}^{(t)}}) \,ln{P(\mathbf{Z}|\mathbf{X},\boldsymbol{\theta}^{(t)})}}_{entropy} \tag{1.4} \end{align} \]
$(1.3)$式構成了$L(\boldsymbol{\theta})$的下界，而$(1.4)$式的右邊爲$P(\mathbf{Z}|\mathbf{X},\boldsymbol{\theta}^{(t)})的熵 \geqslant 0$ ，其獨立於咱們想要優化的參數$\boldsymbol{\theta}$，於是是一個常數。因此極大化$L(\boldsymbol{\theta})$的下界$(1.3)$式就等價於極大化$\mathcal{Q}(\boldsymbol{\theta}, \boldsymbol{\theta}^{(t)})$，$\mathcal{Q}(\boldsymbol{\theta}, \boldsymbol{\theta}^{(t)})$ (Q函數) 亦可表示爲 $\,\mathbb{E}_{\mathbf{Z}|\mathbf{X},\boldsymbol{\theta}^{(t)}}\,lnP(\mathbf{X},\mathbf{Z}|\boldsymbol{\theta})$，其完整定義以下：

基於觀測數據 $\mathbf{X}$ 和當前參數$\theta^{(t)}$計算未觀測數據 $\mathbf{Z}$ 的條件機率分佈$P(\mathbf{Z}|\mathbf{X}, \theta^{(t)})$，則Q函數爲徹底數據的對數似然函數關於$\mathbf{Z}$的指望。

此即E步中指望值的來歷。

接下來來看M步。在$(1.3)$式中若令$\boldsymbol{\theta} = \boldsymbol{\theta}^{(t)}$，則下界$(1.3)$式變爲：
\[ \begin{align*} & \sum\limits_{\mathbf{Z}}P(\mathbf{Z|\mathbf{X},\boldsymbol{\theta}^{(t)}}) \,ln\frac{P(\mathbf{X},\mathbf{Z}|\boldsymbol{\theta}^{(t)}) }{P(\mathbf{Z}|\mathbf{X},\boldsymbol{\theta}^{(t)})} \\ =\;\; & \sum\limits_{\mathbf{Z}}P(\mathbf{Z|\mathbf{X},\boldsymbol{\theta}^{(t)}}) \,ln\frac{P(\mathbf{Z|\mathbf{X},\boldsymbol{\theta}^{(t)}})P(\mathbf{X}|\boldsymbol{\theta}^{(t)})}{P(\mathbf{Z}|\mathbf{X},\boldsymbol{\theta}^{(t)})} \\ = \;\; & \sum\limits_{\mathbf{Z}}P(\mathbf{Z|\mathbf{X},\boldsymbol{\theta}^{(t)}}) \,lnP(\mathbf{X}|\boldsymbol{\theta}^{(t)}) \\ = \;\; & lnP(\mathbf{X}|\boldsymbol{\theta}^{(t)}) \;\;=\;\; L(\boldsymbol{\theta}^{(t)}) \end{align*} \]
能夠看到在第t步，$L(\boldsymbol{\theta}^{(t)})$的下界與$L(\boldsymbol{\theta}^{(t)})$相等，又因爲極大化下界與極大化Q函數等價，於是在M步選擇一個新的$\boldsymbol{\theta}$來極大化$\mathcal{Q}(\boldsymbol{\theta}, \boldsymbol{\theta}^{(t)})$，就能使$L(\boldsymbol{\theta}) \geqslant \mathcal{Q}(\boldsymbol{\theta}, \boldsymbol{\theta}^{(t)}) \geqslant \mathcal{Q}(\boldsymbol{\theta}^{(t)}, \boldsymbol{\theta}^{(t)}) = L(\boldsymbol{\theta}^{(t)})$ (這裏爲了便於理解就將$\mathcal{Q}(\boldsymbol{\theta}, \boldsymbol{\theta}^{(t)})$與$(1.3)$式等同了)，也就是說$L(\boldsymbol{\theta})$是單調遞增的，經過EM算法的不斷迭代能保證收斂到局部最大值。

EM算法流程：

輸入：觀測數據$\mathbf{X}$，隱變量$\mathbf{Z}$，聯合機率分佈$P(\mathbf{X},\mathbf{Z}|\boldsymbol{\theta})$

輸出：模型參數$\boldsymbol{\theta}$

初始化參數$\boldsymbol{\theta}^{(0)}$
E步：利用當前參數$\boldsymbol{\theta}^{(t)}$計算Q函數， $\mathcal{Q}(\boldsymbol{\theta}, \boldsymbol{\theta}^{(t)}) = \sum\limits_{\mathbf{Z}}P(\mathbf{Z|\mathbf{X},\boldsymbol{\theta}^{(t)}}) \,ln{P(\mathbf{X},\mathbf{Z}|\boldsymbol{\theta})}$
M步：極大化Q函數，求出相應的 $\boldsymbol{\theta} = \mathop{argmax}\limits_{\boldsymbol{\theta}}\mathcal{Q}(\boldsymbol{\theta}, \boldsymbol{\theta}^{(t)})$
重複 2. 和3. 步直至收斂。

EM算法也可用於極大後驗估計，極大後驗估計僅僅是在極大似然估計的基礎上加上參數$\boldsymbol{\theta}$的先驗分佈，即 $p(\boldsymbol{\theta})p(\mathbf{X}|\boldsymbol{\theta})$，則取對數後變爲$ln\,p(\mathbf{X}|\boldsymbol{\theta}) + ln\,p(\boldsymbol{\theta})$，因爲後面的$ln\,p(\boldsymbol{\theta})$不包含隱變量$\mathbf{Z}$，因此E步中求Q函數的步驟不變。而在M步中須要求新的參數$\mathbf{\theta}$，所以須要包含這一項，因此M步變爲
\[ \boldsymbol{\theta} = \mathop{argmax}\limits_{\boldsymbol{\theta}} \left[\mathcal{Q}(\boldsymbol{\theta}, \boldsymbol{\theta}^{(t)}) + ln(p(\boldsymbol{\theta})\right] \]

$\large{\S} \normalsize\mathrm{3}$ EM算法深刻

上一節中遺留了一個問題：爲何式$(1.2)$中引入的分佈是$P(\mathbf{Z}|\mathbf{X},\boldsymbol{\theta}^{(t)})$而不是其餘分佈？下面以另外一個角度來闡述。

假設一個關於隱變量$\mathbf{Z}$的任意分佈$q(\mathbf{Z})$，則運用指望值的定義，$(1.1)$式變爲：
\[ \begin{align*} L(\boldsymbol{\theta}) = lnP(\mathbf{X}|\boldsymbol{\theta}) &= \sum\limits_{\mathbf{Z}}q(\mathbf{Z})\,lnP(\mathbf{X}|\boldsymbol{\theta}) \quad\qquad \text{上下同乘以 $q(\mathbf{Z}) \,P(\mathbf{X},\mathbf{Z}|\boldsymbol{\theta})$}\\ & = \sum\limits_{\mathbf{Z}}q(\mathbf{Z}) ln\frac{P(\mathbf{X}|\boldsymbol{\theta})q(\mathbf{Z}) P(\mathbf{X},\mathbf{Z}|\boldsymbol{\theta})}{q(\mathbf{Z}) P(\mathbf{X},\mathbf{Z}|\boldsymbol{\theta})} \\ & = \sum\limits_{\mathbf{Z}}q(\mathbf{Z}) ln\frac{P(\mathbf{X},\mathbf{Z}|\boldsymbol{\theta})}{q(\mathbf{Z})} + \sum\limits_{\mathbf{Z}}q(\mathbf{Z}) ln \frac{P(\mathbf{X}|\boldsymbol{\theta})q(\mathbf{Z}) }{P(\mathbf{X},\mathbf{Z}|\boldsymbol{\theta})} \\ & = \sum\limits_{\mathbf{Z}}q(\mathbf{Z}) ln\frac{P(\mathbf{X},\mathbf{Z}|\boldsymbol{\theta})}{q(\mathbf{Z})} + \sum\limits_{\mathbf{Z}}q(\mathbf{Z}) ln \frac{q(\mathbf{Z}) }{P(\mathbf{Z}|\mathbf{X},\boldsymbol{\theta})} \\ & = \underbrace{\sum\limits_{\mathbf{Z}}q(\mathbf{Z}) ln\frac{P(\mathbf{X},\mathbf{Z}|\boldsymbol{\theta})}{q(\mathbf{Z})}}_{L(q,\boldsymbol{\theta})} + \mathbb{KL}(q(\mathbf{Z})||P(\mathbf{Z}|\mathbf{X},\boldsymbol{\theta}))) \tag{2.1} \end{align*} \]
$(2.1)$式的右端爲$q(\mathbf{Z})$和後驗分佈$P(\mathbf{Z}|\mathbf{X},\boldsymbol{\theta})$的KL散度，由此 $lnP(\mathbf{X}|\boldsymbol{\theta})$被分解爲$L(q,\boldsymbol{\theta})$和$\mathbb{KL}(q||p)$ 。因爲KL散度總大於等於0，因此$L(q,\boldsymbol{\theta})$是$lnP(\mathbf{X}|\boldsymbol{\theta})$的下界，如圖：

由此可將EM算法視爲一個座標提高(coordinate ascent)的方法，分別在E步和M步不斷提高下界$L(q,\boldsymbol{\theta})$，進而提高$lnP(\mathbf{X}|\boldsymbol{\theta})$ 。

在E步中，固定參數$\boldsymbol{\theta}^{old}$，當且僅當$\mathbb{KL}(q||p) = 0$，即$L(q,\boldsymbol{\theta}) = lnP(\mathbf{X}|\boldsymbol{\theta})$時，$L(q,\boldsymbol{\theta})$達到最大，而$\mathbb{KL}(q||p) = 0$的條件是$q(\mathbf{Z}) = P(\mathbf{Z}|\mathbf{X}, \boldsymbol{\theta})$，所以這就是式$(1.2)$中選擇分佈$P(\mathbf{Z}|\mathbf{X},\boldsymbol{\theta}^{old})$的緣由，如此一來$L(q,\boldsymbol{\theta})$ 也就與$(1.3)$式一致了。

在M步中，固定分佈$P(\mathbf{Z}|\mathbf{X},\boldsymbol{\theta}^{old})$，選擇新的$\boldsymbol{\theta}^{new}$來極大化$L(q,\boldsymbol{\theta})$ 。同時因爲$P(\mathbf{Z}|\mathbf{X},\boldsymbol{\theta}^{old}) \neq P(\mathbf{Z}|\mathbf{X},\boldsymbol{\theta}^{new})$，因此$\mathbb{KL}(P(\mathbf{Z}|\mathbf{X},\boldsymbol{\theta}^{old}) || P(\mathbf{Z}|\mathbf{X},\boldsymbol{\theta}^{new})) > 0$，致使$lnP(\mathbf{X}|\boldsymbol{\theta})$提高的幅度會大於$L(q,\boldsymbol{\theta})$提高的幅度，如圖：

所以在EM算法的迭代過程當中，經過交替固定$\boldsymbol{\theta}$ 和 $P(\mathbf{Z}|\mathbf{X},\boldsymbol{\theta}^{old})$來提高下界$L(q,\boldsymbol{\theta})$ ，進而提高對數似然函數$L(\boldsymbol{\theta})$ ，從而在隱變量存在的狀況下實現了極大似然估計。在下一篇中將探討EM算法的具體應用。

Reference: