Armadillo 线性代数库中的聚类算法避坑

    1、本文的由来

         最近由于需要在C++语言编写的项目中使用高斯混合模型聚类算法，最开始是打算自己写一个的（参考的是《机器学习》，周志华著这本书），但是最后发现自己写的算法运行效率低，而且对于维度比较高的样本集（样本维度超过5），会出现聚类失败的情况。原因是高维度的聚类计算过程中会出现浮点溢出的问题，计算过程浮点溢出是高斯混合模型的天然缺陷。因为公式中存在大量的乘法和指数运算，这些运算会使得数据被迅速放大或者缩小，很快就会达到现在计算机提供的数据类型所能表达的数据极限，从而造成数据溢出。下面给出高斯模型的公式：

 这是单个高斯模型公式，高斯混合模型则是多个这样的公式相加在一起形成，相加的时候每个高斯模型前面还会有一个权重系数，高斯混合模型公式如下：

 表示高斯模型的数量，在聚类中表示分类数，表示第g个高斯模型的权重系数。有的朋友可能会问为啥相加时前面要加一个权重系数呢，对于这个问题我也纳闷过，这个大概要搞数学的人才能给出完美的解析。但是如果对非线性曲线拟合或者线性曲线拟合有了解的话，其实你会发现高斯混合模型算法就是高级复杂一些的非线性拟合，和其他拟合一样都是尝试通过已知的样本，去预估混合模型中的各个参数的最优值。只是这里求最优解的方程变成了似然函数（对数似然函数），用的方法也不是最小二乘法、梯度下降这类的，而是EM算法。和其他拟合一样，当你给一个样本集选择一个拟合模型（函数模型）的时候，也就默认了样本集合是可以被这个函数模型所表达出来的。这里也一样，当你选择了高斯混合模型来聚类一个样本集的时候，你就默认了高斯混合模型是可以正确表达你的样本集合的，但是问题就出在这里了，我们怎么确定对于给定的样本集合高斯混合模型是可以正确对他进行表达的呢，没有人能告诉你答案，而且高斯混合模型大部分时候也不能100%正确的对给定的样本集合进行表达，特别是样本量很大且维度很高的时候，我们能做的只是让高斯混合模型尽可能准确地表达给定的样本集合。也就是说高斯混合模型的应用过程中会出现各种各样的不确定性（模型是否适合样本集？），那么给每个高斯模型加上一个权重系数对不确定因素进行表达，拟合出来的模型或许会更符合样本集的实际情况。以上是本人对权重系数的理解，可能不正确。如果没看明白这一段也没有关系，这对后续的阅读没有影响。

        以上是高斯混合模型的基础公式的讲解，这部分内容网上很多地方都有很好的讲解，这里仅从整体的角度说一下，更多的细节可以网上找别人的博客，这里不再多说。个人建议如果你要在你的项目中用到高斯混合模型聚类，最好还是要对算法原理有一定的了解，不然，很有可能在使用过程中会出现聚类效果不理想，你又不知道怎么解决的情况，本人就是活生生的例子。

2、Armadillo 线性代数库中的高斯混合模型聚类算法

         自己实现算法不太行之后，我便找第三方的聚类算法库。网上找了一圈，也没找到有几个库支持高斯混合模型聚类算法，主要是很少有人会用C++语言做聚类算法分析，大都选用Python或者java等其他高级语言实现了，用C++做聚类分析相关功能的估计不多。Armadillo库的下载安装之类的这里不再多说，网上很多教程，这里给出Armadillo库的学习文档网址：Armadillo: C++ library for linear algebra & scientific computing 要用好这个库，参考官网的学习文档是必不可少的。

        Armadillo 线性代数库中与高斯混合模型聚类算法相关的有两个类：gmm_diag、gmm_full。这两个类都可以做高斯混合模型聚类，主要的区别在于gmm_diag用的协方差矩阵是对角线矩阵（对角线上的元素不为0,其他元素为0），gmm_full用的协方差矩阵是全矩阵，对应完整的高斯模型公式。gmm_diag可以理解为gmm_full的化简版本。那么什么时候该用gmm_diag什么时候该用gmm_full呢？当样本的各个维度值是相互独立的，不会相互影响的时候，那么用gmm_diag是比较好的选择，不仅可以减少计算量，而且可以避免模型出现病态解，一般情况使用这个类就可以了。当样本的各个维度值不是相互独立的，会相互影响的时候，那么用gmm_full是比较好的选择。关于这两个类使用的算法，Armadillo库的学习文档里面有一个链接可以查看，不过是英文版的，链接如下：http://arma.sourceforge.net/armadillo_spcs_2017.pdf

        但是无论你选择使用gmm_diag和gmm_full中的哪一个，都有可能会出现病态解，病态解的出现是求解高斯混合模型时所用的算法导致的，几乎无法消除，目前几乎所有的聚类算法或者无监督类的学习算法都只能求得局部最优解，而局部最优有的时候是病态的解。你想求全局最优解？那你要自己发明一个可以做到这一点的算法了，目前好像没有办法可以做到。所以算法解出来的模型参数很有可能和我们想要的聚类模型参数相去甚远，这就会导致聚类结果不准确。即便是用k-means聚类算法也会出现聚类结果不准确（不是我们想要的结果）的情况，我刚开始用Armadillo 库中gmm_full::learn()方法的时候就出现了病态解，那是对一个7维的样本集做聚类，聚类的结果简直垃圾，还不如直接用k-means算法，一开始是以为调用learn函数没给对参数或者样本点质量不好导致的聚类不准确，后来才发现是高斯混合模型得到了病态解，那么如何解决这个问题呢？想要彻底消除病态解那是不可能的，但是我们可以通过一些现有的办法去避免病态解出现的可能。就针对Armadillo库中的gmm_diag和gmm_full类而言，有效避免出现病态解的方法有两个（我会的只有这两个）：

2.1 两个避免病态解的方法

1）手动设定初始质心，不要使用随机初始质心。首先无论是高斯混合聚类算法还是k-means聚类算法，在算法的开始的时候，都需要选择一个初始的质心。有了初始质心，算法才能根据已有质心开始一轮迭代，得到新的质心，再用新的质心继续迭代求解下一个质心，如此重复，直到满足算法结束条件。所以初始质心的选定，会直接影响算法的聚类结果。如果随机选择初始质心，那么聚类结果也具有一定的随机性，时好时坏。要避免这种情况，可以使用其他方法选出初始质心，初始质心之间的差别要尽量大一些，即：质心间的区分度要好。这里可以使用最大最小距离聚类算法，初步确定质心的位置，最大最小距离聚类算法可以参考这个链接：聚类算法-最大最小距离算法(实例+代码)_pan_jinquan的博客-CSDN博客_最大最小距离聚类算法

使用最大最小距离聚类算法求出质心后，将其设为gmm_diag或gmm_full类的初始质心。代码如下：

arma::gmm_full  gmmFull;//高斯混合模型聚类对象
int             iKmIterNum = 10;//K均值初始化参数时的最大迭代次数
double          dMinCovs = 1e-30;//协方差矩阵对角线上的元素的最小值
int             iDim = 7;//样本的维度
int             iK = 4;//分簇数量
arma::mat       matCentroid;//存储预设的初始质心，可以用最大最小距离聚类算法生成，也可以用其他方法生成，每一列代表一个质心
arma::mat       matSmp;//待聚类样本集合，为实际要用聚类的样本，每一列代表一个样本


m_gmmFull.reset(iDim, iK);//预设特征维度以及分簇数量
m_gmmFull.set_means(matCentroid);//预设聚类的初始均值（初始质心）
bool bRes = m_gmmFull.learn(matSmp, iK, arma::eucl_dist, arma::keep_existing, iKmIterNum, iIterNum, dMinCovs, true);//调用learn方法，完成高斯混合模型的拟合，这里第四个参数一定要选arma::keep_existing这个枚举

2）可以的话尽量使用gmm_diag类做聚类（特别是样本维度比较高的时候，如：超过5维）。因为相较于gmm_full类，gmm_diag类不容易出现病态解。样本维度较低时可以尝试使用gmm_full类，毕竟gmm_full类用的才是完整的高斯混合模型，这样或许可以获得更准确的结果。此外，无论是用这两个类中的哪一个，gmm_xxx.learn()函数中的第7个参数，即上面代码中的dMinCovs参数一定要设置一个比较小的且大于0的值，这个值是协方差矩阵的下限值，算法在运行的过程中会保证协方差矩阵的元素的值不低于这个下限值。有的时候，改变这个值可以直接改变算法的聚类效果。如果遇到病态解，可以尝试改变一下这个参数，或许可以得到你想要的结果。

2.2、如何获取聚类结果以及聚类结果的置信度（准确分类的概率）

相较于k均值算法高斯混合模型聚类算法的另外一个特点是“软分类”，即分簇的结果不是一个非此即彼的分类标签，而是一个描述分簇置信度（可能性）的向量（向量的维度和分簇的数量相同），最终将样本划分给置信度值最大的分簇。那么在Armadillo线性代数库提供的gmm_diag和gmm_full类中如何获取某个样本对应每个分簇的置信度信息呢？非常遗憾，Armadillo库中没有提供直接获取分簇结果置信度相关的函数，只提供了一个gmm_xxx.assign()函数用于获取样本的分簇标签信息以及gmm_xxx.log_p()函数用于获取样本的对数似然概率，但是log_p()函数返回的值不是一个[0,1]之间的概率数，而是一个不确定大小的数，至于为什么不能直接返回概率值，是因为Armadillo库中的高斯混合模型算法是改进过的，用的是对数高斯混合模型，不是原始的高斯混合模型，这么做可以避免算法求解的过程中出现浮点运算溢出，但是也导致了无法直接获得我们想要的置信度（概率）信息。这也是高斯混合模型聚类算法尴尬的地方，本来选择使用高斯混合模型聚类算法就是冲着软分类能提供置信度这个重要的信息来的，结果发现没有办法直接获取这个信息，这和直接用k均值聚类就没有什么区别了，还不如k均值算法效果稳定呢。真是难受！但是理论上来说是可以通过gmm_diag和gmm_full类求得的模型参数算出给定样本数据在各个分簇中的可信度的。以下是我阅读了算法文档后总结出来的求分簇置信度的计算方法，仅供学习参考，这个计算方法不一定正确，毕竟是我自己看官方的算法文档后自己推导出来的，不能保证准确。

方法一：

int Assign_Probability(arma::vec vecSmp, double& dProb)//参数vecSmp是一个列向量，表示一个样本；参数dProb用于输出分簇的置信度（概率值）
{
	int iFlag = m_gmmFull.assign(vecSmp, arma::prob_dist);//获取样本的分簇下标（分簇标签）

	double dAllProb_log = m_gmmFull.log_p(vecSmp);//获取样本的总体对数似然	
	double dProb_log = m_gmmFull.log_p(vecSmp, iFlag) + log(m_gmmFull.hefts(iFlag));//获取样本在分簇下标为iFlag簇的对数似然
	dProb = exp(dProb_log - dAllProb_log);//根据Armadillo库中对数优化后的高斯混合模型公式，反向计算概率值

	return iFlag;
}

方法二：

arma::rowvec Probabilitys(arma::colvec vecSmp)//通过gmm_xxx类拟合出来的高斯混合模型的参数计算样本在各个分簇的置信度；参数vecSmp是一个列向量，表示一个样本。
{
	arma::rowvec  rowProb(m_iK);//用于存储每个高斯模型分量的概率
	arma::mat matX(vecSmp);//将样本向量转成矩阵

	double tmpP = sqrt(pow(2 * arma::datum::pi, m_iDim));//m_iDim为样本维度
	for (register int i = 0;i < m_iK;i++)//m_iK为分簇总数
	{
		double dW = m_gmmFull.hefts(i);//获取当前高斯模型的权重
		arma::mat matMean(m_gmmFull.means.col(i));//获取当前高斯模型的均值向量

		arma::mat matCovs(m_iDim, m_iDim);
		
		matCovs = m_gmmFull.fcovs.slice(i);//对于gmm_full类型
		
		//matCovs.diag() = m_gmmDiag.dcovs.col(i);//对于gmm_diag类型

		double tmpD = sqrt(arma::det(matCovs));
		double tmpH = 1.0 / (tmpP * tmpD);
		arma::mat matRelDist = matX - matMean;//样本和均值的相对距离
		arma::mat matRelDistT = matRelDist.t();//相对距离转置
		arma::mat matTmp = matRelDistT *  matCovs.i() * matRelDist;

		double tmpExp = exp(-0.5 * matTmp(0, 0));
		rowProb(i) = dW * tmpH * tmpExp;
	}
	return rowProb;
}

上述函数中，“m_”前缀开头的变量是我自己写的类的成员变量，这里其实粘贴的是类成员函数的代码，但是我没有重新去定义“m_”前缀开头的变量，而是在代码中注释了这些变量的含义，如果读者需要用到上面的代码，请按照变量的含义自己给出变量的定义。

3、结语

        关于Armadillo 线性代数库中的聚类算法内容就分享到这里，能用到这部分知识的人应该很少。毕竟我也是脑子抽风才想用C++实现聚类功能，可以的话还是用python或者其他高级语言来做聚类的需求才是明智之举。正因为用的人不多，一旦遇到问题就特别蛋疼了（因为没有太多的资料可以借鉴），写这篇文章也是希望能帮助到一些遇到同样问题的同学或朋友。以上内容都是基于本人对高斯混合模型聚类算法的理解来写的，难免有理解不正确或者描述错误的地方，请大家选择性吸收，如发现错误，欢迎指出。最后，创作不易，请读者不要自私转发，转发前应获得作者允许，如未能及时获得作者允许转发的，应注明出处，感谢大家。