Jensen不等式

Jensen不等式給出了積分的凸函數值必定大于凸函數(convex)的積分值的定理。在凸函數曲線上的任意兩點間連接一條線段，那么線段會位于曲線之上，這就是將Jensen不等式應用到兩個點的情況，如圖(1)所示\((t\in[0,1])\)。我們從概率論的角度來描述Jensen不等式：假設\(f(x)\)為關于隨機變量\(x\)的凸函數\(f'(x)\geq 0\)，則有\(f\left(E(x)\right)\leq E\left(f(x)\right)\)。反之，如果\(f(x)\)為關于\(x\)的凹函數(concave)，則有\(f\left(E(x)\right)\geq E\left(f(x)\right)\)。

?

對于嚴格的凸函數\(\left(f''(x)>0\right)\)，想取得相等關系\(f\left(E(x)\right)=E\left(f(x)\right)\)，則有\(x=E(x)\)，那么\(x\)必須為常數。

期望最大化

期望最大化(Expectation Maximization,EM)算法常用于求解概率模型的極大似然問題，其典型的應用場景包括：1)觀察到的數據不完全;2)似然函數無法直接計算但可以用隱變量(latent variables)來表示。對于第一種應用場景，我們可以將缺失的數據視為隱變量，然后用已有的觀測數據將其推導出來，最后也就回到了第二種應用場景。假設在概率模型中，\(x^{(i)}\)和\(z^{(i)}\)分別為觀測值和隱變量(\(i=1,2,\cdots,m\))，兩者的聯合概率為\(p(x,z;\theta)\)，其中\(\theta\)為模型參數。我們的目標是在僅知道\(x^{(i)}\)的前提下，求最佳的模型參數\(\theta\)以最大化似然函數 \begin{equation} \mathcal{L}(\theta)=\sum_{i=1}^m\log\left(x^{(i)};\theta\right)=\sum_{i=1}^m\log\sum_{j=1}^Kp\left(x^{(i)},z^{(i)}=j;\theta\right) \end{equation} 這里簡單假設\(z^{(i)}\)為離散值\(\{1,2,\cdots,K\}\)，在\(z^{(i)}\)為連續值或離散值和連續值的組合時情況也是類似的。 在僅知道\(x^{(i)}\)的情況下用極大似然法估計參數很困難。如果\(z^{(i)}\)只可能取少數幾個離散值，我們還可以根據\(z^{(i)}\)分情況優化似然函數，最終求得最合理的模型參數；但若\(z^{(i)}\)的取值情況很復雜時，前面的方法很顯然是行不通的。既然直接求解似然函數的最大值很難，只能退而求其次，求似然函數下界的最大值，得到一個相對較好的模型參數，這就是EM算法的精髓所在。這里我們引入隱變量\(z^{(i)}\)的概率分布\(Q_i(z^{(i)})\)，則對數似然函數表述為如下形式： \begin{equation} \begin{array}{rl} \mathcal{L}(\theta)&=\sum_{i=1}^m\log\sum_{z^{(i)}=1}^KQ_i\left(z^{(i)}\right)\frac{p\left(x^{(i)},z^{(i)};\theta\right)}{Q_i\left(z^{(i)}\right)}\\ &=\sum_{i=1}^m\log E_{z^{(i)}\sim Q_i}\left[\frac{p\left(x^{(i)},z^{(i)};\theta\right)}{Q_i\left(z^{(i)}\right)}\right] \end{array} \end{equation} \(Q_i(z^{(i)})\)為隱變量\(z^{(i)}\)的概率分布，必須滿足約束條件： \begin{equation} Q_i(z^{(i)})\geq 0,\sum_{z^{(i)}=1}^KQ_i(z^{(i)})=1 \end{equation} 對數函數\(\log x\)在區間\((0,+\infty)\)上為嚴格的凹函數，根據Jensen不等式可得到似然函數的下界： \begin{equation} \begin{array}{rl} \mathcal{L}(\theta)&=\sum_{i=1}^m\log E_{z^{(i)}\sim Q_i}\left[\frac{p\left(x^{(i)},z^{(i)};\theta\right)}{Q_i\left(z^{(i)}\right)}\right]\\ &\geq \sum_{i=1}^mE_{z^{(i)}\sim Q_i}\left[\log \frac{p\left(x^{(i)},z^{(i)};\theta\right)}{Q_i\left(z^{(i)}\right)}\right]=L(\theta) \end{array} \end{equation}

從EM算法的描述可知，EM算法的每次迭代分為兩步：E-step和M-step，如圖(2)所示。

E-step使似然函數下界\(L(\theta)\)在\(\theta^{(t)}\)處緊貼著似然函數\(\mathcal{L}(\theta)\)，我們需要調整\(Q_i(z^{(i)})\)以滿足： \begin{equation} \frac{p\left(x^{(i)},z^{(i)};\theta\right)}{Q_i\left(z^{(i)}\right)}=C \end{equation} 結合\(Q_i(z^{(i)})\)的約束條件，很容易得到如下關系式： \begin{equation} Q_i(z^{(i)})=\frac{p\left(x^{(i)},z^{(i)};\theta\right)}{\sum_{z^{(i)}=1}^K p\left(x^{(i)},z^{(i)};\theta\right)}=\frac{p\left(x^{(i)},z^{(i)};\theta\right)}{p\left(x^{(i)};\theta\right)}=p\left(z^{(i)}|x^{(i)};\theta\right) \end{equation} M-step則求使\(L(\theta^{(t)})\)取得極大值的參數\(\theta^{(t+1)}\)作為下一次迭代的初始點。

如果將\(L(\theta)\)定義為關于\(Q_i(z^{(i)})\)和\(\theta\)的函數\(J(Q,\theta)\)，則EM算法的優化策略和坐標上升(coordinate ascent)是一致的。在E-step中，沿著\(Q\)定義的方向最大化\(J(Q,\theta)\)；在M-step中，沿著\(\theta\)定義的方向最大化\(J(Q,\theta)\)，這也可以解釋EM算法在迭代過程中會使目標函數收斂到局部最優解。

高斯混合模型

密度估計(Density Estimation)中的高斯混合模型（Gaussian Mixture Model,GMM)是EM算法的一個典型應用，下面就以GMM為例展示EM的強大威力。高斯混合模型可以表示成\(K\)個高斯分布的線性加權： \begin{equation} P(x)=\sum_{k=1}^K\phi_k\mathcal{N}(x|\mu_k,\Sigma_k) \end{equation} 其中\(w_k\)為第\(k\)個高斯分布的權值。我們要估算出模型中的所有參數\(\phi_k,\mu_k,\Sigma_k\)。 假設現有觀測數據集合\(\mathcal{S}=\{x^{(i)}|x^{(i)}\in\mathbb{R}^n,i=1,\cdots,m\}\)，每個樣本\(x^{(i)}\)都有一個對應的隱變量\(z^{(i)}\in\{1,\cdots,K\}\)表示\(x^{(i)}\)屬于哪一個高斯分布,\(z^{(i)}\)服從參數為\(\phi_1,\cdots,\phi_K\)的多項分布\(P(z^{(i)}=k)=\phi_k\)，參數滿足以下約束條件： \begin{equation} \sum_{k=1}^K\phi_k=1,\quad \phi_k\geq 0,k=1,\cdots,K \end{equation} \(x^{(i)}\)和\(z^{(i)}\)之間的聯合分布\(P(x^{(i)},z^{(i)})=P(x^{(i)}|z^{(i)})P(z^{(i)})\)。在已知\(z^{(i)}=k\)的前提下我們知道\(x^{(i)}\)的概率分布，即\(x^{(i)}|z^{(i)}\sim \mathcal{N}(\mu_k,\Sigma_k)\): \begin{equation} P(x^{(i)}|z^{(i)}=k)=\frac{1}{(2\pi)^{\frac{n}{2}}|\Sigma_k|^{\frac{1}{2}}}\exp\left(-\frac{1}{2}(x^{(i)}-\mu_k)^T\Sigma_k^{-1}(x^{(i)}-\mu_k)\right) \end{equation} 很可惜\(z^{(i)}\)是未知的隱變量，可能是任意一個高斯分布，因此\(P(x^{(i)})=\sum_{z^{(i)}=1}^K P(x^{(i})|z^{(i)})P(z^{(i)})\)。似然函數形式如下： \begin{equation} \ell({\phi,\mu,\Sigma})=\prod_{i=1}^mP(x^{(i)};\phi,\mu,\Sigma)=\prod_{i=1}^m\left(\sum_{z^{(i)}=1}^K P(x^{(i}|z^{(i)})P(z^{(i)})\right) \end{equation} 將上式轉化為形式上更簡單的對數似然函數并加上約束條件構成拉格郎日函數： \begin{equation} \mathcal{L}({\phi,\mu,\Sigma})=\sum_{i=1}^m\log\left(\sum_{z^{(i)}=1}^K P(x^{(i}|z^{(i)})P(z^{(i)})\right)+\lambda(\sum_{k=1}^K\phi_k-1) \end{equation} 仔細觀察這個拉格郎日函數，我們發現對數函數的輸入為多個高斯函數的加權，這就導致了利用普通的求偏導方法是無法獲得一個解析解的。根據前面的EM算法，我們先執行E-step，利用貝葉斯定律猜測每個樣本服從每個高斯分布的概率： \begin{equation} Q_k^{(i)}=P(z^{(i)}=k|x^{(i)})=\frac{P(x^{(i)}|z^{(i)}=k)P(z^{(i)}=k)}{\sum_{k=1}^KP(x^{(i)}|z^{(i)}=k)P(z^{(i)}=k)} \end{equation} 然后求得與上述拉格郎日函數等價的下界函數\(L(\phi,\mu,\Sigma)\): \begin{equation} \begin{array}{rl} L(\phi,\mu,\Sigma)&=\sum_{i=1}^m\sum_{k=1}^KQ_k^{(i)}\log P(x^{(i)},z^{(i)}=k)+\lambda(\sum_{k=1}^K\phi_k-1)\\ &=\sum_{i=1}^m\sum_{k=1}^KQ_k^{(i)}\left(\log P(x^{(i)}|z^{(i)}=k)+\log P(z^{(i)}=k)\right)\\ &\quad+\lambda(\sum_{k=1}^K\phi_k-1) \end{array} \end{equation} 注意到該下界函數中的對數函數已經直接作用到每個高斯分布上了。現在我們可以對參數求偏導來計算使下界函數最大的更好的參數。 \begin{equation} \frac{\partial L}{\partial\mu_k}=0\Rightarrow\mu_k=\frac{\sum_{i=1}^mQ_k^{(i)}x^{(i)}}{\sum_{i=1}^mQ_k^{(i)}} \end{equation} \begin{equation} \frac{\partial L}{\partial\Sigma_k}=0\Rightarrow\Sigma_k=\frac{\sum_{i=1}^mQ_k^{(i)}(x^{(i)}-\mu_k)(x^{(i)}-\mu_k)^T}{\sum_{i=1}^mQ_k^{(i)}} \end{equation} \begin{equation} \frac{\partial L}{\partial\phi_k}=0\Rightarrow\sum_{i=1}^m\frac{Q_k^{(i)}}{\phi_k}+\lambda=0 \end{equation} 結合約束條件\(\sum_{k=1}^K\phi_k=1\)，得到\(\lambda=-m\),因此\(\phi_k=\sum_{i=1}^mQ_k^{(i)}/m\)。 根據上述推到，我們給出GMM模型的算法： ?EM算法收斂速度要比K-means慢，而且每次迭代的計算量也很大。較好的初始值可以加快收斂速度，因此通常可以用K-means完成粗略的聚類，然后用每個類簇的均值、方差和類簇占總樣本的比例來初始化GMM的參數。似然函數可能存在多個局部最優解，且EM算法肯定可以搜索到似然函數的局部最大值，但不能保證可以求得全局的最大值。

根據上述的GMM算法，我寫了一個 Python版本的GMM代碼 。我隨機生成了3組高斯分布的數據，并在圖(3)中用三種不同顏色的統計直方圖表示；然后用生成的數據用EM算法訓練GMM模型，訓練過程中的對數似然函數的變化曲線如圖(4)所示；最后用訓練好的GMM模型反過來預測訓練數據的概率密度函數值，如圖(3)中的紅色曲線所示。

Expectation Maximization and GMM