混合高斯聚类的一种实现方法

近来看了一篇文章《Hierarchical Clustering of a Mixture Model》【1】，它是一篇比较早的文章了，2005年，Jacob Goldberger Sam Roweis，Department of Computer Science, University of Toronto。文中讲到一种简单的高斯聚类方法，近来有文章将它用于目标检测（Object Detection）最后阶段，以代替NMS（non maximum suppression，非极大值抑制），获得更好的bbox预测。我试了一下，这种方法有其独到之处，经它处理后的bbox相互重叠的情况有很大改善，而且较NMS要“软”一些，没那么“硬”，生成的box样子要比NMS能更好地覆盖目标。接下来，我将结合自己的代码和实验进行小结。

一、高斯聚类

混合高斯模型（MoG，Gaussian Mixture Model）是一种常见的参数化概率模型，其表达形式如下：
$$f(y)=\sum _{i=1}^k{\alpha }_{i}N(y;{\mu }_i,{\Sigma }_i)=\sum _{i=1}^k{\alpha }_i{f}_i(y)(1)$$
$f(y)$是由 $k$ 个d维高斯分布构成的混合分布，各高斯分布 $f_i(y)$ 的期望和方差分别为 ${\mu }_i,{\Sigma }_i$。可将每个独立的高斯分布称为一个高斯核，所谓高斯聚类就是将多个高斯核进行聚类，用较少的高斯核来近似表达它，此过程描述如下：
$$f(y)=\sum _{i=1}^k{\alpha }_i{f}_i(y)\approx \sum _{j=1}^m{\beta }_j{g}_j(y)\text{, and }k>m(2)$$
(2)式中 $g_j(y)$ 也是与$f_i(y)$ 同维度的高斯分布，且约等于式左边的高斯核数量要大于右边的高斯核数量。所谓高斯聚类，指的就是用较少的高斯核混合分布来拟合较多核的高斯混合分布。
要完成 $G(y)={\sum }_{j=1}^m{\beta }_j{g}_j(y)$ 对 $F(y)={\sum }_{i=1}^k{\alpha }_i{f}_i(y)$ 的拟合即求使分布 $G(y)$ 与分布 $F(y)$ 距离最小的参数，设 $\theta$ 是 $G(y)$ 的可调参数集，于是 $G(y)$ 可写为 $G_{\theta }(y)$，于是拟合问题就是：
$$\theta =\arg \underset{\theta }{\min }Distance(G_{\theta },F)(3)$$
$G(y)$ 与 $F(y)$ 是两个分布，衡量概率分布的距离可以用KL散度：
$$KL(G,F)={\mathbb{E}}_F(\log \frac{p_g}{p_f})(4)$$
有了衡量拟合效果的距离定义（KL），有了可控参数模型（$G(y)$），似乎拟合问题就可以直接转化为最优化问题了：将 $G_{\theta }(y)$ 和 $F(y)$ 是代入（3）、（4），距离KL对参数集（θ={ β,μ,Σ}\theta=\{\beta,\mu,\Sigma\}θ={ β,μ,Σ}）求偏导，偏导置零求解最优参数集。但不幸的是$G_{\theta }(y)$ 和 $F(y)$ 都是混合高斯模型，这个过程没有闭式解，直接求最优解的方法行不通。怎么办？
突破点1：
两个混合高斯的KL没有闭式，但两个高斯的KL却是有闭式形式的：
设两个高斯分布分别为 $N_1({\mu }_1,{\Sigma }_1)$ 和 ${\mathbb{N}}_2({\mu }_2,{\Sigma }_2)$，则它们的KL距离为：
$$KL(N_1||N_2)=\frac{1}{2}\left(\log \frac{|{\Sigma }_2|}{|{\Sigma }_1|}+Tr({\Sigma }_2^{-1}{\Sigma }_1)+({\mu }_1-{\mu }_2{)}^T({\Sigma }_2{)}^{-1}({\mu }_1-{\mu }_2)+d\right)$$
证明：
d维高斯分布：
N(x∣u,Σ)=1(2π)n/2∣Σ∣1/2exp⁡{ −12(x−u)TΣ−1(x−u)} N(x|u,\Sigma ) = {1 \over { { {(2\pi )}^{n/2}}{ {\left| \Sigma \right|}^{1/2}}}}\exp \{ - {1 \over 2}{(x - u)^T}{\Sigma ^{ - 1}}(x - u)\} N(x∣u,Σ)=(2π)n/2∣Σ∣1/21​exp{ −21​(x−u)TΣ−1(x−u)}
则两个高斯分布的KL散度为：
DKL(P1∣∣P2)=EP1[log⁡P1−log⁡P2] =12EP1[−log⁡∣Σ1∣−(x−u1)TΣ1−1(x−u1)+log⁡∣Σ2∣+(x−u2)TΣ2−1(x−u2)] =12log⁡∣Σ2∣∣Σ1∣+12EP1[−(x−u1)TΣ1−1(x−u1)+(x−u2)TΣ2−1(x−u2)] =12log⁡∣Σ2∣∣Σ1∣+12EP1{ −tr[Σ1−1(x−u1)(x−u1)T]+tr[Σ2−1(x−u2)(x−u2)T]}