自动分层混合PPCA

注意：该博客来自于论文：A Hierarchical Latent Variable Model for Data Visualization和Automated Hierarchical Mixtures of Probabilistic Principal Component Analyzers，仅用于学习，转载请注明出处。

1 混合PPCA

使用有限混合模型聚类是一种众所周知的方法。在本模型中，假定数据是从混合分量的密度函数中生成的，其中每一个分量的混合密度函数$p(\mathbf{t}\mid i)$代表一个cluster。PPCA是一个概率模型，我们用一个single PPCA分布来建立每一个混合的分量。因此观测变量$\mathbf{t}$的概率密度可写成：
$$p(\mathbf{t})=\sum _{i=1}^{k_0}{\pi }_{i}p\left(\mathbf{t}\mid {\mu }_i,{\sigma }_i^2,{\mathbf{W}}_i\right)$$
其中，$p\left(\mathbf{t}\mid {\mu }_i,{\sigma }_i^2,{\mathbf{W}}_i\right)$代表第$i$个component的PPCA密度函数，$k_0$代表component的数量，${\pi }_i$表示第$i$个component的混合比例，所有的${\pi }_i$服从一个约束：${\pi }_i\ge 0$和${\sum }_{i=1}^{k_0}{\pi }_i=1$。

观测数据的对数似然为：

$$\mathcal{L}=\sum _{n=1}^N\log \left\{\sum _{i=1}^{k_0}{\pi }_{i}p\left({\mathbf{t}}_n\mid {\mu }_i,{\sigma }_i^2,{\mathbf{W}}_i\right)\right\}$$

但是，上述的对数似然很难直接优化，因此我们使用EM算法来寻找局部最大值。如果我们假定一组示性变量$z_{ni}$(通常为缺失数据),他们具体指定由哪个模型负责生成哪个数据点${\mathbf{t}}_n$，然后完全数据的对数似然为：

$${\mathcal{L}}_c=\sum _{n=1}^N\sum _{i=1}^{k_0}{z}_{ni}\log \left\{{\pi }_{i}p\left({\mathbf{t}}_n\mid {\mu }_i,{\sigma }_i^2,{\mathbf{W}}_i\right)\right\}$$

应用EM算法来计算参数${\pi }_i,{\mu }_i,{\sigma }_i^2$和${\mathbf{W}}_i$的MLE：

E-step:

属于第$i$个component的数据点${\mathbf{t}}_n$的后验概率$R_{ni}$为：
$$R_{ni}=E\left[z_{ni}\right]=\frac{{\pi }_{i}p\left({\mathbf{t}}_n\mid {\mu }_i,{\sigma }_i^2,{\mathbf{W}}_i\right)}{p\left({\mathbf{t}}_n\right)}$$

M-step:

$$\begin{array}{c}{\tilde{\pi }}_i=\frac{1}{N}{\sum }_{n=1}^N{R}_{ni}\\ {\tilde{\mu }}_i=\frac{{\sum }_{n=1}^N{R}_{ni}{\mathbf{t}}_n}{{\sum }_{n=1}^N{R}_{ni}}\end{array}$$
为了更新${\sigma }_i^2$和${\mathbf{W}}_i$，首先计算加权样本协方差矩阵：
$${\mathbf{S}}_i=\frac{{\sum }_{n=1}^N{R}_{ni}\left({\mathbf{t}}_n-{\tilde{\mu }}_i\right){\left({\mathbf{t}}_n-{\tilde{\mu }}_i\right)}^T}{{\sum }_{n=1}^N{R}_{ni}}$$
然后计算：
$$\begin{array}{c}{\sigma }_{ML}^2=\frac{1}{d-q}{\sum }_{i=q+1}^d{\lambda }_i\\ {\mathbf{W}}_{ML}={\mathbf{U}}_q{\left({\Lambda }_q-{\sigma }_{ML}^2\mathbf{I}\right)}^{1/2}\mathbf{R}\end{array}$$

2 分层混合PPCA

将混合PPCA模型扩展为分层混合PPCA模型。下面的表述也将是非常一般地。我们可以将此分层混合模型应用到任意混合参数密度模型中。

截止目前为止，我们考虑的模型都是：第一层为一个简单的潜变量模型，第二层为$k_0$个component模型的混合。现在将这样的两层混合模型扩展到三层混合模型，那么第三层次一组潜变量模型$G_i$是由第二层次的第$i$个component模型生成的。那么对应的概率密度可以表述为：
$$p(\mathbf{t})=\sum _{i=1}^{k_0}{\pi }_i\sum _{j\in g_i}{\pi }_{j\mid i}p\left(\mathbf{t}\mid {\mu }_{i,j},{\sigma }_{i,j}^2,{\mathbf{W}}_{i,j}\right)$$
其中，$p\left(\mathbf{t}\mid {\mu }_{i,j},{\sigma }_{i,j}^2,{\mathbf{W}}_{i,j}\right)$代表一个single PPCA模型，${\pi }_{j\mid i}$代表的是混合比例，其中，${\pi }_{j\mid i}\ge 0$ and ${\sum }_j{\pi }_{j\mid i}=1$。

每一个层次对应一个生成模型，更低的水平代表更精细更更详细的表示。

注意：图片来自Bishop的论文：A Hierarchical Latent Variable Model
for Data Visualization

第三层模型参数的确定同样可以再次被视为数据缺失问题，其中缺失的信息对应于指定哪个模型生成每个数据点的标签。所以当未提供标签信息时，对数似然为：

$$L=\sum _{n=1}^N\ln \left\{\sum _{i=1}^{k_0}{\pi }_i\sum _{j\in g_i}{\pi }_{j\mid i}p({\mathbf{t}}_n\mid {\mu }_{i,j},{\sigma }_{i,j}^2,{\mathbf{W}}_{i,j})\right\}$$

然而，如果我们给出一组示性变量$z_{ni}$，指定每个数据点$t_n$是由第二层上的第$i$个模型生成的，那么对数似然就会变成：
$$L=\sum _{n=1}^N\sum _{i=1}^{k_0}{z}_{ni}\log \left\{{\pi }_i\sum _{j\in g_i}{\pi }_{j\mid i}p\left({\mathbf{t}}_n\mid {\mu }_{i,j},{\sigma }_{i,j}^2,{\mathbf{W}}_{i,j}\right)\right\}$$

事实上，对于每个模型$i$生成数据点$t_n$的后验概率$R_{ni}$形式，我们只有从层次结构的第2层获得部分的概率的信息。对上述公式求期望，得到第三层次水平的对数似然：

$$L=\sum _{n=1}^N\sum _{i=1}^{k_0}{R}_{ni}\log \left\{{\pi }_i\sum _{j\in g_i}{\pi }_{j\mid i}p\left({\mathbf{t}}_n\mid {\mu }_{i,j},{\sigma }_{i,j}^2,{\mathbf{W}}_{i,j}\right)\right\}$$

上述式子同样可以使用EM算法进行最大化，除了E-STEP中模型$(i,j)$生成数据点${\mathbf{t}}_n$的后验概率之外，这和一般的混合模型有相同的形式。
$$R_{ni,j}=R_{ni}{R}_{nj\mid i}$$
其中，
$$R_{nj\mid i}=\frac{{\pi }_{j\mid i}p\left({\mathbf{t}}_n\mid {\mu }_{i,j},{\sigma }_{i,j}^2,{\mathbf{W}}_{i,j}\right)}{{\sum }_{k\in g_i}{\pi }_{k\mid i}p\left(\mathbf{t}\mid {\mu }_{i,k},{\sigma }_{i,k}^2,{\mathbf{W}}_{i,k}\right)}$$
M-step中，混合系数和均值的更新公式为：

$$\begin{array}{l}{\tilde{\pi }}_{j\mid i}=\frac{{\sum }_n{R}_{ni,j}}{{\Sigma }_n{R}_{ni}}\\ {\tilde{\mu }}_{i,j}=\frac{{\sum }_n{R}_{ni,j}{\mathbf{t}}_n}{{\sum }_n{R}_{mi,j}}\end{array}$$

缺失数据的$z_{ni,j}$的后验期望为：
$$\begin{array}{c}⟨{\mathbf{x}}_{ni,j}⟩={\mathbf{M}}_{i,j}^{-1}{\mathbf{W}}_{i,j}^{\mathrm{T}}\left({\mathbf{t}}_n-{\tilde{\mu }}_{i,j}\right)\\ ⟨{\mathbf{x}}_{ni,j}{\mathbf{x}}_{ni,j}^{\mathrm{T}}⟩={\sigma }_{i,j}^2{\mathbf{M}}_{i,j}^{-1}+⟨{\mathbf{x}}_{ni,j}⟩⟨{\mathbf{x}}_{ni,j}^{\mathrm{T}}⟩\end{array}$$
$$\begin{array}{c}{\tilde{\mathbf{W}}}_{i,j}=[{\sum }_n{R}_{ni,j}\left({\mathbf{t}}_n-{\tilde{\mathbf{\mu }}}_{i,j}\right)⟨{\mathbf{x}}_{ni,j}^{\mathrm{T}}⟩{\left[{\sum }_n{R}_{ni,j}⟨{\mathbf{x}}_{ni,j}{\mathbf{x}}_{ni,j}^{\mathrm{T}}⟩\right]}^{-1}\\ {\tilde{\sigma }}_{i,j}^2=\frac{1}{d{\sum }_n{R}_{ni,j}}\{{\sum }_n{R}_{ni,j}{\Vert {\mathbf{t}}_n-{\tilde{\mu }}_{i,j}\Vert }^2\\ -2{\sum }_n{R}_{ni,j}⟨{\mathbf{x}}_{ni,j}^{\mathrm{T}}\backslash {\tilde{\mathbf{W}}}_{i,j}^{\mathrm{T}}\left({\mathbf{t}}_n-{\tilde{\mathbf{\mu }}}_{i,j}\right)\\ +{\sum }_n{R}_{ni,j}\mathrm{Tr}\left[⟨{\mathbf{x}}_{ni,j}{\mathbf{x}}_{ni,j}^{\mathrm{T}}⟩{\tilde{\mathbf{W}}}_{i,j}^{\mathrm{T}}{\tilde{\mathbf{W}}}_{i,j}\right]\}\end{array}$$
我们可以将E-step替换进M-step方程来获得更新公式：

$$\begin{array}{c}{\tilde{\mathbf{W}}}_{i,j}={\mathbf{S}}_{i,j}{\mathbf{W}}_{i,j}{\left({\sigma }_{i,j}^2\mathbf{I}+{\mathbf{M}}_{i,j}^{-1}{\mathbf{W}}_{i,j}^{\mathrm{T}}{\mathbf{S}}_{i,j}{\mathbf{W}}_{i,j}\right)}^{-1}\\ {\tilde{\sigma }}_{i,j}^2=\frac{1}{d}\mathrm{Tr}\left({\mathbf{S}}_{i,j}-{\mathbf{S}}_{i,j}{\mathbf{W}}_{i,j}{\mathbf{M}}_{i,j}^{-1}{\tilde{\mathbf{W}}}_{i,j}^{\mathrm{T}}\right)\end{array}$$
其中，
$${\mathbf{S}}_{i,j}=\frac{1}{N_{i,j}}\sum _n{R}_{ni,j}\left({\mathbf{t}}_n-{\tilde{\mathbf{\mu }}}_{i,j}\right){\left({\mathbf{t}}_n-{\tilde{\mathbf{\mu }}}_{i,j}\right)}^{\mathrm{T}}$$
$$N_{i,j}=\sum _n{R}_{ni,j}$$

${\mathbf{S}}_{i,j}$被解释为responsibility-weighted的协方差矩阵。

对于这种分层建模的技术，可以直接扩展到其他任何的参数族。

3 自动分层混合PPCA

3.1 建立分层和决定什么时候进行分裂clusters

我们按照如下构造分层结构：
在每一层上，我们对每一个cluster执行一个阶次设别检验来观察是否该cluster应该在下一层分割成两个子clusters。特别地，我们采用分层的混合PPCA模型，$g_i$中有两个子clusters。然后，我们基于一种标准评价的测度来比较原来的cluster和它的两个子clusters，如果两个子clusters的结果优于原来的cluster，那么我们就用这两个子clusters来替换原来的cluster，如果原来的cluster优于两个子clusters，那么原来的cluster就不能分割成两个子clusters，我们直接将原来的cluster直接复制到第三层。我们只考虑三层模型，不考虑更低层次的模型了。我们重复这个过程，直到第三层所有的单个cluster都不能分割成两个子clusters或者clusters的数量达到最大$K_{max}$为止。

这引导我们：应该采用什么样的标准来进行分割呢？？

这个问题类似于如何确定混合模型中的clusters的数量，这是一个比较困难的任务。我们使用积分分类似然(ICL)准则：

$$ICL={\mathcal{L}}_c-m\log (N)/2$$

${\mathcal{L}}_c$是上述提到的完全的对数似然，$m$代表需要估计的模型的参数的个数，$N$代表样本量，$k$代表混合的component的个数，$m$随$k$的变化而变化。

在此，我们选择用$R_{ni}$来替换$z_{ni}$，这代表混合问题中‘soft’聚类解。这个标准成为 ICL-BICin (McLachlan & Peel, 2000) 并且他们演示这个标准优于BIC和AIC准则。

通过最大化ICL来选择clusters的数量。ICL看起来很像BIC，但是它并不是惩罚观测对数似然，而是惩罚完全数据对数似然函数的期望${\mathcal{L}}_c$,${\mathcal{L}}_c$的期望等于$\mathcal{L}+{\sum }_{n=1}^N{\sum }_{i=1}^{k_0}{R}_{ni}\log R_{ni}$(也就是似然函数减去模糊分类矩阵的估计熵)。ICL测量的是观测数据对数似然减去clusters重叠程度再减去模型参数复杂性的惩罚。

让我们总结一下我们用于计算单个PPCA组分及其两个子类ICL的表达式:

第$i$个component的ICL为：

$$\sum _{n=1}^N{R}_{ni}\log \left\{{\pi }_{i}p\left({\mathbf{t}}_n\mid {\mu }_i,{\sigma }_i^2,{\mathbf{W}}_i\right)\right\}-\frac{m_1\log (N)}{2}$$

第$i$个component的两个子类的ICL为：
$$\begin{array}{r}{\sum }_{n=1}^N{\sum }_{j=1}^2{R}_{ni,j}\log \left\{{\pi }_{i,j}p\left({\mathbf{t}}_n\mid {\mu }_{i,j},{\sigma }_{i,j}^2,{\mathbf{W}}_{i,j}\right)\right\}-\frac{m_2\log (N)}{2}\end{array}$$

$m_1$和$m_2$表示第$i$个component及其两个子类的自有参数的个数。

3.2 决定每个component的主子空间的维度

我们介绍了一种类似于传统主成分分析降维技术的简单而快速的方法。

$$q_i={\mathrm{argmin}}_{1<q<d}\left(\frac{{\sum }_{j=1}^q{\lambda }_j}{\mathrm{trace}\left(S_i\right)}>90\%\right)$$
${\lambda }_i$是$S_i$的特征值。