本篇将介绍整理各种聚类算法，包括k-means，GMM(Guassian Mixture Models, 高斯混合)，EM(Expectation Maximization，期望最大法)，Spectral Clustering(谱聚类)，Mean Shift(均值偏移)和DBSCAN(Density-Based Spatial Clustering of Applications with Noise)

0. 先验的基础知识

文章: 聚类算法的先验基础知识

1. K-Means

K-means算法是一种常用的聚类算法，用于将数据集中的样本分成 K 个不同的类别或簇。

• 原理介绍：

  • 初始化： 随机选择 K 个初始聚类中心。
  • 分配样本： 将每个样本分配到距离最近的聚类中心所在的簇。
  • 更新中心： 对每个簇，计算其所有样本的均值，将该均值作为新的聚类中心。
  • 重复步骤2和步骤3，直到收敛： 当聚类中心不再变化或变化很小时，算法收敛，得到最终的聚类结果。

• 优缺点：

  • 优点：
    • 简单直观，易于理解和实现。
    • 可以处理大规模数据集。
    • 对处理球形簇形状的数据效果较好。
  • 缺点：
    • 需要预先指定聚类个数 K。
    • 对异常值敏感，可能导致聚类结果不佳。
    • 结果可能收敛到局部最优解，因此多次运行可以获得更好的结果。
  • 应用：
    • 数据分析：对数据集进行聚类分析，发现数据集中的内在结构。
    • 图像处理：对图像像素进行聚类，实现图像分割和压缩。
    • 推荐系统：对用户进行聚类，提供个性化推荐。

• C++代码：

#include <iostream>
#include <vector>
#include <cmath>
#include <limits>using namespace std;// 计算两个点之间的欧氏距离
double distance(vector<double>& point1, vector<double>& point2) {double sum = 0.0;for (int i = 0; i < point1.size(); ++i) {sum += pow(point1[i] - point2[i], 2);}return sqrt(sum);
}// 更新簇的中心
void updateCenters(vector<vector<double>>& points, vector<vector<double>>& centers, vector<int>& clusterAssignments) {vector<int> clusterCounts(centers.size(), 0);vector<vector<double>> newCenters(centers.size(), vector<double>(centers[0].size(), 0.0));for (int i = 0; i < points.size(); ++i) {int cluster = clusterAssignments[i];for (int j = 0; j < points[i].size(); ++j) {newCenters[cluster][j] += points[i][j];}clusterCounts[cluster]++;}for (int i = 0; i < newCenters.size(); ++i) {if (clusterCounts[i] > 0) {for (int j = 0; j < newCenters[i].size(); ++j) {newCenters[i][j] /= clusterCounts[i];}}}centers = newCenters;
}// K-means算法
void kMeans(vector<vector<double>>& points, int k, double threshold) {vector<vector<double>> centers(k);vector<int> clusterAssignments(points.size(), -1);// 随机初始化聚类中心for (int i = 0; i < k; ++i) {int randIndex = rand() % points.size();centers[i] = points[randIndex];}bool converged = false;while (!converged) {converged = true;// 分配样本到最近的聚类中心for (int i = 0; i < points.size(); ++i) {double minDist = numeric_limits<double>::max();int minIndex = -1;for (int j = 0; j < centers.size(); ++j) {double dist = distance(points[i], centers[j]);if (dist < minDist) {minDist = dist;minIndex = j;}}if (clusterAssignments[i] != minIndex) {converged = false;clusterAssignments[i] = minIndex;}}// 更新聚类中心if (!converged) {vector<vector<double>> oldCenters = centers;updateCenters(points, centers, clusterAssignments);// 判断中心点变化是否小于阈值double maxCenterShift = 0.0;for (int i = 0; i < centers.size(); ++i) {double shift = distance(oldCenters[i], centers[i]);if (shift > maxCenterShift) {maxCenterShift = shift;}}if (maxCenterShift < threshold) {converged = true; // 达到收敛条件}}}// 打印聚类结果for (int i = 0; i < points.size(); ++i) {cout << "Point " << i << " belongs to cluster " << clusterAssignments[i] << endl;}
}int main() {// 示例数据vector<vector<double>> points = {{1, 2}, {5, 6}, {3, 4}, {8, 9}, {7, 8}};int k = 2; // 聚类个数double threshold = 0.001; // 中心点位置变化阈值// 调用 K-means 算法kMeans(points, k, threshold);return 0;
}

• 技巧:
  因为在更新中心点之前需要得到所有点所属的中心点，然后才能计算更新后的中心点。这样效率比较低，可以每次迭代只选择其中的一部分点来完成中心点的更新，这种做法叫做mini batch k_means。

2. GMM

高斯混合模型（Gaussian Mixture Model，GMM）是一种常用的概率模型，用于对数据进行聚类或密度估计。原理是基于假设：数据是由多个高斯分布组成的混合物。每个高斯分布代表一个聚类簇，而混合模型则将多个高斯分布加权组合起来，以适应数据中的不同分布情况。简单来说就是集合P中包含的所有点是来自多个分布的，但是现在想用所有分布的加权平均来表达集合P中所有点的分布。但是因为单一的高斯分布不明确，且加权平均权重也不明确，所以通过这组数据得到$K$个不同的高斯分布以及$K$个权重就是GMM聚类算法的核心。最后的聚类结果就是得到每个点属于每个高斯分布的后验概率，就完成了聚类。

说人话：有一组2维点，表示为$X=[{\mathbf{x}}_1,{\mathbf{x}}_2,...,{\mathbf{x}}_n]，{\mathbf{x}}_i\in R_2$，${\mathbf{x}}_i$表示第$i$个点的坐标。这$n$个点来自于$K$个不同的高斯分布(也就是说这$n$个点可以分成$K$组，每组点都是从同一个分布中采样得到)，但是我们得到的信息只有这组点的坐标，并不知道这K个高斯分布参数是多少。另外，K个单一的高斯分布对每个点的概率贡献权重也不知道。所以GMM的算法核心就是通过这组点的坐标信息得到$\pi ,\mu ,\Sigma ,\pi \in R_{2\times k},\mu \in R_{2\times k},k=1,2,3,...,K$。$\Sigma =[{\Sigma }_1,{\Sigma }_2,...,{\Sigma }_k],{\Sigma }_i=diag({\sigma }_{i1}^2,{\sigma }_{i2}^2);$这里高斯分布的参数之所以都是2维的，是因为我们的点坐标是2维的，这里代表不同维度有不同的高斯分布参数，比如说颜色属性和坐标属性。

单维度变量$x$的高斯分布，描述为:
N ( x ∣ μ , σ 2 ) = 1 ( 2 π σ 2 ) 1 2 e x p { − 1 2 σ 2 ( x − μ ) 2 } \mathcal{N}(x|\mu,\sigma^2) = {1 \over (2\pi\sigma^2)^{1\over2}} exp \{{-1\over2\sigma^2}(x-\mu)^2\} N(x∣μ,σ2)=(2πσ2)21​1​exp{2σ2−1​(x−μ)2}

多维度变量$\mathbf{x},\mathbf{x}\in R_D$的高斯公式，描述为：
$$\mathcal{N}(\mathbf{x}\mid \mathbf{\mu },\mathbf{\Sigma })=\frac{1}{(2\pi {)}^{D/2}}\frac{1}{|\mathbf{\Sigma }{|}^{1/2}}\exp \left\{-\frac{1}{2}(\mathbf{x}-\mathbf{\mu }{)}^{\mathrm{T}}{\mathbf{\Sigma }}^{-1}(\mathbf{x}-\mathbf{\mu })\right\}$$

高斯混合分布，公式表达为:
$$p(x)=\sum _{k==1}^K{\pi }_{k}N(x|{\mu }_k,{\Sigma }_k)$$
公式表示了变量$x$在高混合分布中的概率值。

• 算法步骤

  • 初始化：选择初始均值$\mu$、协方差矩阵$\Sigma$和混合系数$\pi ,\mu =[{\mu }_1,{\mu }_2,...,{\mu }_k],\Sigma =[{\Sigma }_1,{\Sigma }_2,...,{\Sigma }_k],\pi =[{\pi }_1,{\pi }_2,...,{\pi }_k],k=1,2,...,K$。
  • EM算法：迭代进行E步（Expectation）和M步（Maximization）。在E步中，计算每个数据点$\mathbf{x}$属于每个高斯分布的后验概率$P(z_k|\mathbf{x}),z_k$表示第$k$个高斯分布。在M步中，利用这些后验概率更新高斯分布的参数$\mu$、$\Sigma$和$\pi$。
  • 收敛：当参数变化小于设定阈值或达到最大迭代次数时停止迭代。

• E-Step，计算每个数据点$x$属于每个高斯分布的后验概率$P(z_k|\mathbf{x})$
  $$\begin{array}{rl}{\gamma }_{nk}\equiv p\left(z_k\mid {\mathbf{x}}_n\right)& =\frac{{\pi }_{k}p\left({\mathbf{x}}_n\mid z_k\right)}{{\sum }_{j=1}^K{\pi }_{j}p\left({\mathbf{x}}_n\mid z_j\right)}\\ & =\frac{{\pi }_k\mathcal{N}\left({\mathbf{x}}_n\mid {\mu }_k,{\Sigma }_k\right)}{{\sum }_{j=1}^K{\pi }_j\mathcal{N}\left({\mathbf{x}}_n\mid {\mu }_j,{\Sigma }_j\right)}\end{array}$$
  $r_{nk}$表示第n个点属于第k个高斯分布的概率。

• M-Step，使用E步求出的后验概率更新高斯分布的参数$\mu$、$\Sigma$和$\pi$。
  首先，需要明确的一点是，更新是需要等式的，那么$r_{nk}$与三者的联系是什么，就需要进一步的求解推导。这里借助极大似然估计的思想，拉格朗日算子法的手段得到三者的联系。
  核心思想是寻找使得观测到的数据在给定模型下出现的可能性最大的参数值。即，选择使似然度函数$L(\mu ,\Sigma ,\pi |x)$最大化的参数值$\mu ,\Sigma ,\pi$。公式表达为:
  $$\pi ,\mu ,\mathbf{\Sigma }=\underset{\pi ,\mu ,\mathbf{\Sigma }}{\mathrm{argmax}}\ln p(\mathbf{X}\mid \pi ,\mu ,\mathbf{\Sigma })=\underset{\pi ,\mu ,\mathbf{\Sigma }}{\mathrm{argmax}}\sum _{n=1}^N\ln \left\{\sum _{k=1}^K{\pi }_k\mathcal{N}\left({\mathbf{x}}_n\mid {\mu }_k,{\mathbf{\Sigma }}_k\right)\right\},s.t.=\sum _{k=1}^K{\pi }_k=1$$

  • 对${\mu }_k$求偏导得，
    $$\begin{gathered} 0=-\sum _{n==1}^N{\gamma }_{nk}{\Sigma }_k({\mathbf{x}}_n-{\mu }_k) \\ {\mu }_k=\frac{1}{N_k}\sum _{n=1}^N{\gamma }_{nk}{\mathbf{x}}_n \end{gathered}$$
    其中，$N_k$是后验概率之和，物理意义是表示有多少点是属于第k个高斯分布的，${\mu }_k$物理意义是属于第k个高斯分布的所有点的加权平均，不是平均值，还有权重

  • 对${\Sigma }_k$求偏导得，
    $${\Sigma }_k=\frac{1}{N_k}\sum _{n-1}^N\gamma \left(z_{nk}\right)\left({\mathbf{x}}_n-{\mu }_k\right){\left({\mathbf{x}}_n-{\mu }_k\right)}^T$$

  • 对拉格朗日参数$\lambda$求偏导得，这里结合拉格朗日算子法，目标函数为：
    ln ⁡ p ( X ∣ π , μ , Σ ) + λ ( ∑ k = 1 K π k − 1 ) ln ⁡ p ( X ∣ π , μ , Σ ) = ∑ n = 1 N ln ⁡ { π k N ( x n ∣ μ k , Σ k ) } \ln p(\mathbf{X} \mid \pi, \mu, \Sigma)+\lambda\left(\sum_{k=1}^{K} \pi_{k}-1\right) \\ \ln p(\mathbf{X} \mid \pi, \mu, \Sigma)=\sum_{n=1}^{N} \ln \left\{\pi_{k} \mathcal{N}\left(\mathbf{x}_{\mathbf{n}} \mid \mu_{k}, \Sigma_{k}\right)\right\} lnp(X∣π,μ,Σ)+λ(k=1∑K​πk​−1)lnp(X∣π,μ,Σ)=n=1∑N​ln{πk​N(xn​∣μk​,Σk​)}
    得，
    $$\lambda =-N$$

  • 对${\pi }_k$求偏导得，
    $${\pi }_k=\frac{N_k}{N}$$
    其中N为所有点的数量

• 代码部分

#include <iostream>
#include <vector>
#include <Eigen/Dense> // 使用Eigen库进行矩阵运算using namespace std;
using namespace Eigen;// 定义GMM模型参数
struct GMMParameters {vector<double> weights; // 每个高斯分布的权重vector<VectorXd> means; // 每个高斯分布的均值向量vector<MatrixXd> covariances; // 每个高斯分布的协方差矩阵
};// 计算高斯分布的概率密度
double gaussianPDF(VectorXd x, VectorXd mean, MatrixXd covariance) {int dim = x.size();double det = covariance.determinant();MatrixXd invCov = covariance.inverse();VectorXd diff = x - mean;double exponent = -0.5 * diff.transpose() * invCov * diff;double normalizer = pow(2 * M_PI, -dim / 2) * sqrt(det);return exp(exponent) / normalizer;
}// 对每个数据点计算属于每个高斯分布的后验概率
MatrixXd computePosterior(vector<VectorXd>& data, GMMParameters& params) {int numData = data.size();int numComponents = params.weights.size();MatrixXd posterior(numData, numComponents);for (int i = 0; i < numData; ++i) {VectorXd x = data[i];for (int j = 0; j < numComponents; ++j) {double likelihood = gaussianPDF(x, params.means[j], params.covariances[j]);posterior(i, j) = params.weights[j] * likelihood;}posterior.row(i) /= posterior.row(i).sum(); // 归一化后验概率}return posterior;
}// 更新GMM参数
void updateGMMParameters(vector<VectorXd>& data, MatrixXd& posterior, GMMParameters& params) {int numData = data.size();int numComponents = params.weights.size();VectorXd weights(numComponents);vector<VectorXd> means(numComponents);vector<MatrixXd> covariances(numComponents);// 更新每个高斯分布的权重、均值和协方差矩阵for (int j = 0; j < numComponents; ++j) {weights[j] = posterior.col(j).sum() / numData; // 更新权重VectorXd mean = VectorXd::Zero(data[0].size());for (int i = 0; i < numData; ++i) {mean += posterior(i, j) * data[i];}mean /= posterior.col(j).sum(); // 更新均值means[j] = mean;MatrixXd covariance = MatrixXd::Zero(data[0].size(), data[0].size());for (int i = 0; i < numData; ++i) {VectorXd diff = data[i] - mean;covariance += posterior(i, j) * (diff * diff.transpose());}covariance /= posterior.col(j).sum(); // 更新协方差矩阵covariances[j] = covariance;}// 更新参数params.weights = weights;params.means = means;params.covariances = covariances;
}int main() {// 生成示例数据vector<VectorXd> data;data.push_back(VectorXd(2));data[0] << 1, 2;data.push_back(VectorXd(2));data[1] << 3, 4;// 初始化GMM参数GMMParameters params;params.weights = {0.5, 0.5};params.means.push_back(VectorXd(2));params.means[0] << 1, 2;params.means.push_back(VectorXd(2));params.means[1] << 3, 4;params.covariances.push_back(MatrixXd::Identity(2, 2));params.covariances.push_back(MatrixXd::Identity(2, 2));// EM算法迭代更新参数for (int iter = 0; iter < 10; ++iter) {MatrixXd posterior = computePosterior(data, params);updateGMMParameters(data, posterior, params);}// 输出更新后的GMM参数cout << "Updated GMM Parameters:" << endl;cout << "Weights: " << params.weights[0] << ", " << params.weights[1] << endl;cout << "Mean 1: " << params.means[0].transpose() << endl;cout << "Mean 2: " << params.means[1].transpose() << endl;cout << "Covariance 1: " << endl << params.covariances[0] << endl;cout << "Covariance 2: " << endl << params.covariances[1] << endl;return 0;
}

3. EM算法

K-means和GMM算法都属于EM算法，EM算法的通用推导比较麻烦，就不再推导了。
实在是写不动了，先贴图吧，以后有空了编辑

4.Spectral Clustering

Spectral Clustering（谱聚类）是一种基于图论和谱理论的聚类算法，其原理、应用场景以及优缺点如下：

• 原理步骤

  • 构建相似度图： 首先根据数据集中样本点之间的相似度构建一个相似度图（通常是无向加权图），其中节点表示样本点，边的权重表示样本点之间的相似度。
  • 计算拉普拉斯矩阵： 根据相似度图计算拉普拉斯矩阵（Laplacian Matrix）。常用的有度数矩阵（Degree Matrix）和邻接矩阵（Adjacency Matrix）构建拉普拉斯矩阵。
  • 特征值分解： 对拉普拉斯矩阵进行特征值分解，得到特征值和对应的特征向量。
  • 选取特征向量： 根据特征值和特征向量，选取其中特征值较小的前 k 个特征向量（k 为聚类簇的数量）。
  • 聚类： 利用选取的特征向量进行聚类，常用的方法包括 K-means 等聚类算法。

• 应用场景

  • 图像分割： 在图像处理中，谱聚类可以用于图像分割，将图像中相似的像素点聚类到同一个分割区域。
  • 社交网络分析： 在社交网络分析中，可以利用谱聚类对用户进行聚类，发现社区结构或用户群体。
  • 模式识别： 在模式识别领域，谱聚类可用于将相似模式或特征聚类到一起，用于数据降维或分类。
  • 推荐系统： 在推荐系统中，可以利用谱聚类对用户或物品进行聚类，实现个性化推荐。

• 优缺点

  • 优点：

1. 对非线性、非凸分布的数据具有较好的适应性。
2. 能够发现不规则形状的聚类簇。
3. 对噪声数据较为鲁棒。

• 缺点：

4. 对于大规模数据集计算复杂度较高。
5. 需要事先确定聚类簇的数量。
6. 对于高维稀疏数据表现不佳。
7. 总体来说，谱聚类算法在处理非线性、非凸分布的数据，发现不规则形状的聚类簇方面具有一定优势，但在大规模数据和
   高维稀疏数据上可能存在一些挑战。因此，在选择聚类算法时需要根据具体应用场景和数据特点进行综合考虑。

谱聚类关注的是连接性，并不是工作在欧式空间上

5. Mean Shift

一种基于密度的聚类算法

1. 初始化
   从数据集中随机选择一个数据点作为初始的圆心。
   设定初始的圆半径$r$
2. 移动圆
   将圆移动到窗口内的数据点的中心位置。
   计算窗口内所有数据点的平均位置，将圆心移动到该位置。
3. 收敛判断
   检查圆心移动的距离是否小于设定的阈值 $ϵ$，如果是则认为已经收敛，停止迭代。
   否则，继续迭代执行步骤2。
4. 移除重叠圆
   如果存在重叠的圆，则保留包含数据点数量最多的圆，移除其他重叠的圆。
   对于重叠的圆，可以比较它们所包含的数据点数量来决定保留哪一个圆。
5. 确定聚类中心
   使用保留的圆的圆心作为聚类中心。
   可以根据圆的圆心与数据点的距离来确定数据点所属的聚类。
6. 重复步骤 1-5
   继续选择新的圆心并执行移动、收敛判断、移除重叠圆、确定聚类中心的步骤，直到所有圆都收敛且没有重叠。
   这种方法结合了Mean Shift的密度估计和圆的移动过程，通过移除重叠圆和确定聚类中心的方式来处理重叠的情况，并最终确定聚类结果。这样的方法可以更好地处理重叠的情况，并在一定程度上提高聚类的准确性。

6. DBSCAN

DBSCAN（Density-Based Spatial Clustering of Applications with Noise）是一种基于密度的聚类算法，它可以有效地处理具有任意形状和大小的聚类，并且能够识别和排除噪声数据点。

• DBSCAN算法基于以下两个重要概念：

  • 核心点（Core Point）： 在给定的半径 ϵ 内至少包含指定数量 MinPts 个数据点的点被认为是核心点。
  • 密度可达（Density Reachable）： 如果点 P 在点 Q 的 ϵ 邻域内，并且点 Q 是核心点，则点 P 是密度可达于点 Q。

• 详细步骤

  • 初始化： 设置参数 ϵ 和MinPts，以及数据集。
  • 核心点标记： 对于每个数据点，计算其 ϵ 邻域内的数据点数量，如果大于等于 MinPts，则将该点标记为核心点。
  • 密度可达性： 对于每个核心点，通过递归地检查其 ϵ 邻域内的点是否密度可达，将其加入到同一个簇中。
  • 噪声点处理： 将未被标记为核心点且不属于任何簇的点标记为噪声点。
  • 聚类结果： 最终，所有密度可达的点被划分为不同的簇，噪声点则被排除。

• 优缺点

  • 优点：
    • 能够处理任意形状的聚类，对于密度不均匀的数据效果较好。
    • 能够识别和排除噪声点。
    • 不需要预先指定聚类数量。
  • 缺点：
    • 对参数 ϵ 和 MinPts 比较敏感，需要调优。
    • 在高维数据集上可能会失效，因为密度的定义在高维空间中变得模糊。
    • 对于密度差异较大的数据集，可能会产生较多的噪声点或者无法正确识别聚类