统计算法的理论与实例之二 EM算法实例： 形状约束的非参数混合高斯模型

版权声明：原创作品，欢迎转载，但转载请以超链接形式注明文章来源(planckscale.info)、作者信息和本声明，否则将追究法律责任。

最近比较无聊，干脆把这个系列补完。本文的内容是我表妹之前在芝大读书时的一个小题目，作为EM算法的一个比较复杂的实例来演示是不错的。由于我没找到引文，就不引用了。我不是做统计的，有错误欢迎留言。

\subsubsection{解析推导}
这是一个更为复杂的估计问题实例。要估计的分布是样本空间里一组函数图形对应的分布。具体说来，我们有随机变量\{\boldsymbol{X},Y,Z\},其中\boldsymbol{X}\in \mathbf{R}^p,Y\in \mathbf{R},Z \in \{1,...K\},其中Z是混合高斯模型中的类别标签变量,对于不同的Z值,Y和\boldsymbol{X}存在一个不同的潜在函数关系m_Z(x),其中m_Z(x)为凸函数。概率模型表述为

(1)   \begin{eqnarray*} p(\boldsymbol{X},Y,Z)&=&\sum_j\delta_{Zj}\pi_j\mathcal{N}(m_j(\boldsymbol{X}),\sigma^2_j) \end{eqnarray*}

现在假定Z不可观测，我们从一系列观测值\{\boldsymbol{x}_i,y_i\}中估计出模型的参数\boldsymbol{\theta}\equiv \{\boldsymbol{\pi},\boldsymbol{m(x)},\boldsymbol{\Sigma}=\{\sigma_1,...\sigma_{K}\}\}的真实值\boldsymbol{\theta}_*。

接下来我们用EM算法来处理这一问题。设\mathcal{X}\equiv\{\boldsymbol{x}_{1},...\boldsymbol{x}_N\},\mathcal{Y}\equiv \{y_1,...y_N\}为观测样本,\mathcal{Z}\equiv\{z_1,...z_N\}为观测样本对应的隐变量。我们首先计算p(\mathcal{Z}\vert \mathcal{X},\mathcal{Y},\boldsymbol{\theta}^{(t)})

(2)   \begin{eqnarray*} \gamma_{ij}^{(t)}&\equiv&p(z_i=j\vert \boldsymbol{x}_i, y_i, \boldsymbol{\theta}^{(t)})\\ &=&\frac{p(\boldsymbol{x}_i,y_i,z_i=j \vert \boldsymbol{\theta}^{(t)})}{\sum_k p(\boldsymbol{x}_i,y_i,z_i=k \vert \boldsymbol{\theta}^{(t)})}\\ &=&\frac{\pi_j \mathcal{N}(m_j(\boldsymbol{x}_i),\sigma^2_j)\vert_{y_i}}{\sum_k \pi_k \mathcal{N}(m_k(\boldsymbol{x}_i),\sigma^2_k)\vert_{y_i}} \end{eqnarray*}

类似于边缘分布和条件分布的因子分解性质，可以证明对任意满足q(\mathcal{Z})=\prod_i q_i(z_i)的概率分布q(\mathcal{Z}),有

(3)   \begin{eqnarray*} E_{q(\mathcal{Z})}[\log L(\mathcal{X},\mathcal{Z})]=\sum_i E_{q_i(z_i)}[\log L(\boldsymbol{x}_i,z_i)] \end{eqnarray*}

这意味着Q(\boldsymbol{\theta}\vert\boldsymbol{\theta}^{(t)})可以分解成在每个样本上各自均值之总和的形式

(4)   \begin{eqnarray*} Q(\boldsymbol{\theta}\vert\boldsymbol{\theta}^{(t)}) &\equiv& E_{\mathcal{Z}\vert \mathcal{X},\mathcal{Y},\boldsymbol{\theta}^{(t)}}[\log L(\boldsymbol{\theta};\mathcal{X},\mathcal{Y},\mathcal{Z})]\\ &=& \sum_i E_{z_i\vert \boldsymbol{x}_i,y_i,\boldsymbol{\theta}^{(t)}}[\log L(\boldsymbol{\theta};\boldsymbol{x}_i,y_i,z_i)] \end{eqnarray*}

我们把该性质的证明留到后面。

于是

(5)   \begin{eqnarray*} Q(\boldsymbol{\theta}\vert\boldsymbol{\theta}^{(t)}) &=&\sum_{ij} p(z_i=j\vert \boldsymbol{x}_i,y_i,\boldsymbol{\theta}^{(t)}) \log L(\boldsymbol{\theta};\boldsymbol{x}_i,y_i,z_i=j)\\ &=&\sum_{ij} \gamma_{ij}^{(t)} \log L(\boldsymbol{\theta};\boldsymbol{x}_i,y_i,z_i=j)\\ &=&\sum_{ij} \gamma_{ij}^{(t)} [\log \pi_j - \log\sigma_j - \frac{1}{2}\log{2\pi} - \frac{(y_i-m_j(\boldsymbol{x}_i))^2}{2\sigma_j^2}] \end{eqnarray*}

需要注意的是\boldsymbol{m}(\boldsymbol{x})=\{\boldsymbol{m}_1(\boldsymbol{x}),...\boldsymbol{m}_{K}(\boldsymbol{x})\}包含了无穷个参数需要优化，为了约化这一问题，我们将其离散化，取m_{ji}=m(\boldsymbol{x}_i)为优化参数，仍然简记为\boldsymbol{m}.这样，\boldsymbol{m}(\boldsymbol{x})被约化到用K\times N个参数来表示。

接下来我们需要求解的极大化问题为

(6)   \begin{eqnarray*} \max_{\boldsymbol{\theta},\boldsymbol{\beta}} &&Q(\boldsymbol{\theta}\vert\boldsymbol{\theta}^{(t)})\\ \rm{subject\ to} &&\sum_j \pi_j=1\\ &&m_{ji}-m_{jk} \geq \boldsymbol{\beta}_{jk}^{T} (\boldsymbol{x}_i-\boldsymbol{x}_k)\qquad \forall i,j,k \end{eqnarray*}

这里我们为m(x)加入了凸函数的形状约束。约束的几何意义很明确，即\textbf{函数总是位于每一点的切线之上方}。

由于Q(\boldsymbol{\theta}\vert\boldsymbol{\theta}^{(t)})中\boldsymbol{\pi}的极大化和其他参数相互独立，因而可将问题解耦。问题简化为两个子问题，其中\boldsymbol{\pi}参数对应的子问题为

(7)   \begin{equation*} \left\{ \begin{aligned} \boldsymbol{\pi}^{(t+1)}&=argmax_{\boldsymbol{\pi}} \sum_{ij}\gamma_{ij}^{(t)}\log\pi_j\\ \sum_j \pi_j&=1 \end{aligned} \right. \end{equation*}

可以求出解析解，即

(8)   \begin{eqnarray*} \pi_j^{(t+1)}&=&\frac{\sum_{i}\gamma_{ij}^{(t)}}{\sum_{ij} \gamma_{ij}^{(t)}}=\frac{\sum_i \gamma_{ij}^{(t)}}{N} \end{eqnarray*}

\boldsymbol{m}参数对应的子问题

(9)   \begin{equation*} \left\{ \begin{aligned} \boldsymbol{m}^{(t+1)}&=argmax_{\boldsymbol{m},\boldsymbol{\beta}} \sum_{ij}\gamma_{ij}^{(t)}[- \log\sigma_j - \frac{(y_i-m_{ji})^2}{2\sigma_j^2}]\\ &=\sum_j argmax_{\boldsymbol{m},\boldsymbol{\beta}}[-N\log\sigma_j-\frac{\sum_i \gamma_{ij}^{(t)}(y_i-m_{ji})^2}{2\sigma_j^2}]\\ &m_{jk}-m_{ji} \geq \boldsymbol{\beta}_{ji}^{T} (\boldsymbol{x}_k-\boldsymbol{x}_i)\qquad \forall i,j,k \end{aligned} \right. \end{equation*}

是个二次规划问题。由于\sigma_j,m_{ji}只出现在具有相同j值的项中，这意味着按不同的j值把项分组，每一组都是独立变化的，因而可以作为一个独立的子问题来最大化。可以发现，\textbf{这里的Gauss混合模型中每个component都对应于一个带约束的最小二乘问题}。这一性质可以用来降低计算复杂度，提高算法效率。

在一些情况下，二次规划中的约束可以大幅度的约减。第一种情况是当X为单变量时。不失一般性的，我们假设观测样本已经按照X的大小做升序排列，此时可以用二阶导数非负来作为凸性约束。于是约束可约化为

(10)   \begin{eqnarray*} m_{jk}-m_{ji} &\geq& \beta_{ji} (x_{k}-x_i)\qquad \forall i,j,k\in\{i-1,i+1\}\\ \beta_{j,i+1}-\beta_{j,i}&\geq&0 \qquad\forall i\in\{1...N-1\},j \end{eqnarray*}

这里\beta_{ji}代表了第j个component在第x_i处导数的估计值。需要注意的是，为了得到一点上导数（切线）的估计，我们需要两个约束，约束在左右两侧的函数值都要位于切线上方。

另一种情况则是当X为多维的向量，且有m(\boldsymbol{X})=\sum_\alpha m^\alpha(X_\alpha)时。此时每个分量m^\alpha(X_\alpha)的凸性意味着m(\boldsymbol{X})也是凸的，因此约束可写为

(11)   \begin{eqnarray*} m^{\alpha}_{jk}-m^{\alpha}_{ji} &\geq& \beta^{\alpha}_{ji} (x_{k}-x_i)\qquad \forall \alpha,i,j,k\in\{i-1,i+1\}\\ \beta^{\alpha}_{j,i+1}-\beta^{\alpha}_{j,i}&\geq&0 \qquad\forall \alpha,i\in\{1...N-1\},j \end{eqnarray*}

即对向量函数的每个分量，都满足\eqref{mConstraint1D}所示的约束。

最后我们来证明E_{q(\mathcal{Z})}[\log L(\mathcal{X},\mathcal{Z})]的分解性质，其中q(\mathcal{Z})为任意满足q(\mathcal{Z})=\prod_i q_i(z_i)的概率分布

(12)   \begin{eqnarray*} &&E_{q(\mathcal{Z})}[\log L(\mathcal{X},\mathcal{Z})]\\ &=&\sum_{\mathcal{Z}} q(\mathcal{Z})\log p(\mathcal{X},\mathcal{Z})\\ &=&\sum_i \sum_{\mathcal{Z}} q(\mathcal{Z}) \log p(\boldsymbol{x}_i,z_i) \end{eqnarray*}

可见对每个sample，都需要在分布q(\mathcal{Z})上求\log p(\boldsymbol{x}_i,z_i)的均值。但由于\mathcal{Z}中包含大量无关的z_k\neq z_i,它们都应该被加和为1,从而只剩下q_i(z_i):

(13)   \begin{eqnarray*} &=&\sum_i \sum_{\{z_1,...,z_N\}} \prod_j q_j(z_j) \log p(\boldsymbol{x}_i,z_i)\\ &=&\sum_i \sum_{z_i}\{\underbrace{\sum_{\{z_1,.\slashed z_i.,z_N\}} \prod_{j\neq i} q_j(z_j)}_{=1}\} q_i(z_i) \log p(\boldsymbol{x}_i,z_i)\\ &=&\sum_i \sum_{z_i} q_i(z_i)\log p(\boldsymbol{x}_i,z_i)\\ &=&\sum_i E_{q_i(z_i)}[\log L(\boldsymbol{x}_i,z_i)] \end{eqnarray*}

\subsubsection{算法实现}
算法实现主要需要注意的是防止数值下溢带来的除零问题。技术上，避免概率连乘积，避免指数，让计算的中间结果尽量在对数域中表示是一些基本技巧。

首先，在计算\gamma_{ij}时，可能会出现某个样本点远离所有component的情况，此时分子分母都容易下溢为零，而出现\frac{0}{0}的情况。我们做如下处理

(14)   \begin{eqnarray*} \gamma_{ij}^{(t)}&=&\frac{\pi_j \mathcal{N}(m_j(\boldsymbol{x}_i),\sigma^2_j)\vert_{y_i}}{\sum_k \pi_k \mathcal{N}(m_k(\boldsymbol{x}_i),\sigma^2_k)\vert_{y_i}}\\ &=&\frac{1}{\sum_k \frac{\pi_k}{\pi_j} \frac{\mathcal{N}(m_k(\boldsymbol{x}_i),\sigma^2_k)\vert_{y_i}}{\mathcal{N}(m_j(\boldsymbol{x}_i),\sigma^2_j)\vert_{y_i}}}\\ &=&\frac{1}{\sum_k \frac{\pi_k}{\pi_j} \frac{\sigma_j}{\sigma_k} \exp(-\frac{(y_i-m_k(\boldsymbol{x}_i))^2}{2\sigma_k^2}+\frac{(y_i-m_j(\boldsymbol{x}_i))^2}{2\sigma_j^2})}\\ \end{eqnarray*}

算法中需要计算对数似然函数的上升率，来判断算法的终止条件是否满足。对数似然函数中有大量的概率连乘积，我们作如下处理来消除之

(15)   \begin{eqnarray*} \log L(\boldsymbol{\theta};\mathcal{X},\mathcal{Y}) &=&\log p(\mathcal{X},\mathcal{Y}\vert\boldsymbol{\theta})\\ &=&\log\prod_i p(\boldsymbol{x}_i,y_i\vert\boldsymbol{\theta})\\ &=&\sum_i\log\sum_{z_i} p(\boldsymbol{x}_i,y_i,z_i\vert\boldsymbol{\theta})\\ &=&\sum_i\log \sum_j \frac{\pi_j}{\sigma_j \sqrt{2\pi}} \exp (-\frac{(y_i-m_j(\boldsymbol{x}_i))^2}{2\sigma_j^2})\\ &=&\sum_i\log \frac{\pi_k}{\sigma_k \sqrt{2\pi}} \exp (-\frac{(y_i-m_k(\boldsymbol{x}_i))^2}{2\sigma_k^2})\\ && \sum_j \frac{\pi_j}{\pi_k}\frac{\sigma_k}{\sigma_j} \exp (\frac{(y_i-m_k(\boldsymbol{x}_i))^2}{2\sigma_k^2}-\frac{(y_i-m_j(\boldsymbol{x}_i))^2}{2\sigma_j^2})\\ &=&\sum_i\{\log \frac{\pi_k}{\sigma_k \sqrt{2\pi}} -\frac{(y_i-m_k(\boldsymbol{x}_i))^2}{2\sigma_k^2}\\ && +\log\sum_j \frac{\pi_j}{\pi_k}\frac{\sigma_k}{\sigma_j} \exp (\frac{(y_i-m_k(\boldsymbol{x}_i))^2}{2\sigma_k^2}-\frac{(y_i-m_j(\boldsymbol{x}_i))^2}{2\sigma_j^2})\} \end{eqnarray*}

对于二次规划子问题，我们需要将其化为标准形式，从而可以利用现成的二次规划库来求解。