聚类

聚类是用于获得关于数据结构的直觉的最常见的机器学习中探索数据分析技术之一。它可以被定义为识别数据中的子组的任务，使得相同子组（cluster）中的数据点非常相似，而不同clusters中的数据点非常不同。换句话说，我们试图在数据中找到同质子组，使得每个聚类中的数据点根据相似性度量（例如基于欧几里得距离或基于相关性距离）尽可能相似。使用哪个相似度度量是由特定于应用程序决定的。

与监督学习不同，聚类被认为是一种无监督学习方法，因为我们没有ground truth来将聚类算法的输出与真实标签进行比较以评估其性能。我们只想通过将数据点分组到不同的子组来尝试研究数据的结构。

在这篇文章中，我们将仅介绍Kmeans，由于其简单性，它被认为是最常用的聚类算法之一。

Kmeans算法

Kmeans算法是一种迭代算法，它试图将数据集划分为 Kpre-defined distinct non-overlapping subgroups（clusters），其中每个数据点只属于一个组。它尝试使cluster内数据点尽可能相似，同时使clusters尽可能保持不同（远）。它将数据点分配给一个cluster，使得数据点之间的平方距离与簇的质心（属于该cluster的所有数据点的算术平均值）之和最小。我们在clusters中的变化越小，数据点在同一cluster中的同质性（相似）就越多。

kmeans算法的工作方式如下：

1. 指定聚类数ķ。
2. 通过首先对数据集进行shuffling， 然后在不替换的情况下随机选择K个数据点。
3. 继续迭代，直到质心没有变化。即，向聚类分配数据点不会改变。

• 计算数据点与所有质心之间的平方距离之和。
• 将每个数据点分配给最近的聚类（质心）。
• 通过获取属于每个聚类的所有数据点的平均值来计算聚类的质心。

kmeans所采用的解决问题的方法被称为期望最大化。E-step将数据点分配给最近的聚类。M-step是计算每个聚类的质心。下面是我们如何用数学方法解决它的分解。

目标函数为:

其中wik = 1，表示数据点xi属于簇k ; 否则，wik = 0。此外，μk是xi的cluster的质心。

这是两个部分的最小化问题。我们首先最小化J w.r.t. wik并将μk固定。然后我们最小化J w.r.t. μk并处理wik固定。从技术上讲，我们首先区分J w.r.t. wik并更新聚类分配（E-step）。然后我们区分J w.r.t. μk并在从前一步骤（M-step）进行聚类分配后重新计算质心。因此，E步骤是：

换句话说，将数据点xi分配给最近的cluster，该cluster通过其与cluster的质心的平方距离之和来判断。

M-step是：

这意味着重新计算每个聚类的质心以反映新的分配。

这里有几点需要注意：

由于包括kmeans的聚类算法使用基于距离的测量来确定数据点之间的相似性，因此建议将数据标准化为平均值为零且标准偏差为1，因为几乎总是任何数据集中的特征都具有不同的测量单位例如年龄与收入。

给定kmeans迭代性质和算法开始时质心的随机初始化，不同的初始化可能导致不同的clusters，因为kmeans算法可能陷入局部最优并且可能不会收敛到全局最优。因此，建议使用不同的质心初始化运行算法，并选择产生较低平方距离总和的运行结果。

分配示例不会改变与聚类内编号没有变化是一回事：

实现

我们将在这里使用简单的kmeans实现来说明一些概念。然后我们将使用sklearn更高效的实现为我们处理许多事情。

import numpy as np
from numpy.linalg import norm
 
 
class Kmeans:
 '''Implementing Kmeans algorithm.'''
 
 def __init__(self, n_clusters, max_iter=100, random_state=123):
 self.n_clusters = n_clusters
 self.max_iter = max_iter
 self.random_state = random_state
 
 def initializ_centroids(self, X):
 np.random.RandomState(self.random_state)
 random_idx = np.random.permutation(X.shape[0])
 centroids = X[random_idx[:self.n_clusters]]
 return centroids
 
 def compute_centroids(self, X, labels):
 centroids = np.zeros((self.n_clusters, X.shape[1]))
 for k in range(self.n_clusters):
 centroids[k, :] = np.mean(X[labels == k, :], axis=0)
 return centroids
 
 def compute_distance(self, X, centroids):
 distance = np.zeros((X.shape[0], self.n_clusters))
 for k in range(self.n_clusters):
 row_norm = norm(X - centroids[k, :], axis=1)
 distance[:, k] = np.square(row_norm)
 return distance
 
 def find_closest_cluster(self, distance):
 return np.argmin(distance, axis=1)
 
 def compute_sse(self, X, labels, centroids):
 distance = np.zeros(X.shape[0])
 for k in range(self.n_clusters):
 distance[labels == k] = norm(X[labels == k] - centroids[k], axis=1)
 return np.sum(np.square(distance))
 
 def fit(self, X):
 self.centroids = self.initializ_centroids(X)
 for i in range(self.max_iter):
 old_centroids = self.centroids
 distance = self.compute_distance(X, old_centroids)
 self.labels = self.find_closest_cluster(distance)
 self.centroids = self.compute_centroids(X, self.labels)
 if np.all(old_centroids == self.centroids):
 break
 self.error = self.compute_sse(X, self.labels, self.centroids)
 
 def predict(self, X):
 distance = self.compute_distance(X, old_centroids)
 return self.find_closest_cluster(distance)

应用

kmeans算法非常受欢迎，可用于各种应用，如市场细分，文档聚类，图像分割和图像压缩等。通常我们进行聚类分析的目标是：

• 对我们正在处理的数据结构有一个有意义的直觉。
• 如果我们认为不同子组的行为存在很大差异，则聚类然后预测将为不同子组建立不同模型的位置。这方面的一个例子是将患者聚集到不同的亚组中，并为每个亚组建立模型，以预测心脏病发作风险的概率。

在这篇文章中，我们将在两种情况下应用聚类：

• 间歇泉喷发分割（2D数据集）。
• 图像压缩。

关于间歇泉爆发分割的Kmeans

我们将首先在2D数据集(http://www.kankanyun.com/data/old_faithful.csv)上实现kmeans算法，看看它是如何工作的。该数据集有272个观测值和2个特征。这些数据涵盖了美国怀俄明州黄石国家公园的Old Faithful间歇泉喷发和喷发持续时间之间的等待时间。我们将尝试在数据点中找到K个子组并相应地对它们进行分组。以下是特征说明：

• eruptions （float）：以分钟为单位的喷发时间。
• waiting（int）：等待下一次喷发的时间。

让我们先绘制数据，Python代码如下：

# Modules
import matplotlib.pyplot as plt
from matplotlib.image import imread
import pandas as pd
import seaborn as sns
from sklearn.datasets.samples_generator import (make_blobs,
 make_circles,
 make_moons)
from sklearn.cluster import KMeans, SpectralClustering
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import silhouette_samples, silhouette_score
 
%matplotlib inline
sns.set_context('notebook')
plt.style.use('fivethirtyeight')
from warnings import filterwarnings
filterwarnings('ignore')
 
# Import the data
df = pd.read_csv('../data/old_faithful.csv')
 
# Plot the data
plt.figure(figsize=(6, 6))
plt.scatter(df.iloc[:, 0], df.iloc[:, 1])
plt.xlabel('Eruption time in mins')
plt.ylabel('Waiting time to next eruption')
plt.title('Visualization of raw data');

我们将使用这些数据，因为它是一个二维数据集，所以很容易绘制和可视化地发现聚类。很明显我们有两个clusters。让我们首先对数据进行标准化，并对K=2的标准化数据运行kmeans算法。

# Standardize the data
X_std = StandardScaler().fit_transform(df)
 
# Run local implementation of kmeans
km = Kmeans(n_clusters=2, max_iter=100)
km.fit(X_std)
centroids = km.centroids
 
# Plot the clustered data
fig, ax = plt.subplots(figsize=(6, 6))
plt.scatter(X_std[km.labels == 0, 0], X_std[km.labels == 0, 1],
 c='green', label='cluster 1')
plt.scatter(X_std[km.labels == 1, 0], X_std[km.labels == 1, 1],
 c='blue', label='cluster 2')
plt.scatter(centroids[:, 0], centroids[:, 1], marker='*', s=300,
 c='r', label='centroid')
plt.legend()
plt.xlim([-2, 2])
plt.ylim([-2, 2])
plt.xlabel('Eruption time in mins')
plt.ylabel('Waiting time to next eruption')
plt.title('Visualization of clustered data', fontweight='bold')
ax.set_aspect('equal');

上面的图显示了由它们所属的cluster着色的数据的散点图。在这个例子中，我们选择K=2。符号“*”是每个cluster的质心。我们可以把这两个clusters看作间歇泉在不同的场景下有不同的行为。

接下来，我们将展示中心的不同初始化可能产生不同的结果。我将使用9种不同的random_state来更改中心的初始化并绘制结果。每个绘图的标题将是每个初始化距离的平方之和。

作为补充说明，这个数据集被认为非常简单，并且在不到10次迭代中收敛。因此，为了了解随机初始化对收敛性的影响，我将使用3次迭代来说明这个概念。然而，在实际应用程序中，数据集并不是那么干净和漂亮!

n_iter = 9
fig, ax = plt.subplots(3, 3, figsize=(16, 16))
ax = np.ravel(ax)
centers = []
for i in range(n_iter):
 # Run local implementation of kmeans
 km = Kmeans(n_clusters=2,
 max_iter=3,
 random_state=np.random.randint(0, 1000, size=1))
 km.fit(X_std)
 centroids = km.centroids
 centers.append(centroids)
 ax[i].scatter(X_std[km.labels == 0, 0], X_std[km.labels == 0, 1],
 c='green', label='cluster 1')
 ax[i].scatter(X_std[km.labels == 1, 0], X_std[km.labels == 1, 1],
 c='blue', label='cluster 2')
 ax[i].scatter(centroids[:, 0], centroids[:, 1],
 c='r', marker='*', s=300, label='centroid')
 ax[i].set_xlim([-2, 2])
 ax[i].set_ylim([-2, 2])
 ax[i].legend(loc='lower right')
 ax[i].set_title(f'{km.error:.4f}')
 ax[i].set_aspect('equal')
plt.tight_layout();

如上图所示，我们最终只采用了两种不同的基于不同初始化的聚类方式。我们会选择平方距离最小的那个。

关于图像压缩的Kmeans

在这部分中，我们将实现kmeans来压缩图像。我们将要处理的图像是396 x 396 x 3.因此，对于每个像素位置，我们将有3个8位整数，用于指定红色，绿色和蓝色强度值。我们的目标是将颜色数量减少到30，并仅使用这30种颜色来表示（压缩）照片。要选择要使用的颜色，我们将在图像上使用kmeans算法，并将每个像素视为数据点。这意味着将图像从高度x宽度x通道reshape为（高度*宽度）x通道，即，我们将在三维空间中具有396 x 396 = 156,816个数据点，这是RGB的强度。这样做将允许我们使用每个像素的30个质心来表示图像，并且将图像的尺寸显着减小6倍。原始图像尺寸为396 x 396 x 24 = 3,763，584位; 但是，新的压缩图像将是30 x 24 + 396 x 396 x 4 = 627,984位。巨大的差异来自于我们将使用质心作为像素颜色的查找，并且将每个像素位置的大小减小到4位而不是8位。

从现在开始，我们将使用sklearnkmeans的实现。这里有几点需要注意：

• n_init是使用不同质心初始化运行kmeans的次数。将报告最佳结果。
• tol 是用于声明收敛的cluster内变体度量。
• 默认值init是k-means ++，它可以产生比质心随机初始化更好的结果。

Python代码如下：

# Read the image
img = imread('images/my_image.jpg')
img_size = img.shape
 
# Reshape it to be 2-dimension
X = img.reshape(img_size[0] * img_size[1], img_size[2])
 
# Run the Kmeans algorithm
km = KMeans(n_clusters=30)
km.fit(X)
 
# Use the centroids to compress the image
X_compressed = km.cluster_centers_[km.labels_]
X_compressed = np.clip(X_compressed.astype('uint8'), 0, 255)
 
# Reshape X_recovered to have the same dimension as the original image 128 * 128 * 3
X_compressed = X_compressed.reshape(img_size[0], img_size[1], img_size[2])
 
# Plot the original and the compressed image next to each other
fig, ax = plt.subplots(1, 2, figsize = (12, 8))
ax[0].imshow(img)
ax[0].set_title('Original Image')
ax[1].imshow(X_compressed)
ax[1].set_title('Compressed Image with 30 colors')
for ax in fig.axes:
 ax.axis('off')
plt.tight_layout();

我们可以看到原始图像和压缩图像之间的比较。压缩的图像看起来接近原始图像，这意味着我们能够保留原始图像的大部分特征。对于较少数量的clusters，我们将以较高的压缩率而牺牲图像质量。这种图像压缩方法称为有损数据压缩，因为我们无法从压缩图像重建原始图像。

评估方法

有监督机器学习时，我们可以用ground truth来评估机器学习模型的性能，我们可以用它来评估不同聚类算法的结果，聚类分析没有可靠的评估指标。此外，由于kmeans需要k作为输入并且不从数据中学习它，因此在任何问题中我们应该具有的clusters的数量没有正确的答案。有时领域知识和直觉可能会有所帮助，但通常情况并非如此。在聚类预测方法中，我们可以根据不同的K聚类评估模型的执行情况，因为聚类用于下游建模。

在这篇文章中，我们将介绍两个可能让我们对k有直觉的指标：

• Elbow method
• Silhouette analysis

Elbow method

肘部法给我们一个想法, 一个好的 k 数目的集群将基于平方距离 (SSE) 之间的总和的数据点和他们分配的聚类的质心。我们选择 k 的地方, SSE 开始平坦, 并形成一个elbow。我们将使用间歇泉数据集, 并评估 SSE 的不同值 k, 并观察曲线可能形成elbow和变平的位置。Python代码如下：

# Run the Kmeans algorithm and get the index of data points clusters
sse = []
list_k = list(range(1, 10))
 
for k in list_k:
 km = KMeans(n_clusters=k)
 km.fit(X_std)
 sse.append(km.inertia_)
 
# Plot sse against k
plt.figure(figsize=(6, 6))
plt.plot(list_k, sse, '-o')
plt.xlabel(r'Number of clusters *k*')
plt.ylabel('Sum of squared distance');

上图显示k=2不是一个好的选择。有时仍然很难计算出大量的簇，因为曲线是单调递减的，可能没有任何elbow，或者有一个明显的点，曲线开始变平。

Silhouette analysis

Silhouette analysis可以用来确定簇之间的分离程度。对每个样本:

• 计算同一cluster(ai)中所有数据点的平均距离。
• 计算距离最近的cluster(bi)中所有数据点的平均距离。
• 计算系数（系数可以取区间[-1,1]中的值）:

• 如果它是0 - >样本非常接近相邻的cluster。
• 它是1 - >样本远离相邻的cluster。
• 它是-1 - >样本被分配给错误的cluster。

因此，我们希望系数尽可能大并且接近1以具有良好的cluster。我们将再次使用间歇泉数据集，因为运行Silhouette analysis的成本更低，实际上很可能只有两组数据点。Python代码如下：

for i, k in enumerate([2, 3, 4]):
 fig, (ax1, ax2) = plt.subplots(1, 2)
 fig.set_size_inches(18, 7)
 
 # Run the Kmeans algorithm
 km = KMeans(n_clusters=k)
 labels = km.fit_predict(X_std)
 centroids = km.cluster_centers_
 
 # Get silhouette samples
 silhouette_vals = silhouette_samples(X_std, labels)
 
 # Silhouette plot
 y_ticks = []
 y_lower, y_upper = 0, 0
 for i, cluster in enumerate(np.unique(labels)):
 cluster_silhouette_vals = silhouette_vals[labels == cluster]
 cluster_silhouette_vals.sort()
 y_upper += len(cluster_silhouette_vals)
 ax1.barh(range(y_lower, y_upper), cluster_silhouette_vals, edgecolor='none', height=1)
 ax1.text(-0.03, (y_lower + y_upper) / 2, str(i + 1))
 y_lower += len(cluster_silhouette_vals)
 
 # Get the average silhouette score and plot it
 avg_score = np.mean(silhouette_vals)
 ax1.axvline(avg_score, linestyle='--', linewidth=2, color='green')
 ax1.set_yticks([])
 ax1.set_xlim([-0.1, 1])
 ax1.set_xlabel('Silhouette coefficient values')
 ax1.set_ylabel('Cluster labels')
 ax1.set_title('Silhouette plot for the various clusters', y=1.02);
 
 # Scatter plot of data colored with labels
 ax2.scatter(X_std[:, 0], X_std[:, 1], c=labels)
 ax2.scatter(centroids[:, 0], centroids[:, 1], marker='*', c='r', s=250)
 ax2.set_xlim([-2, 2])
 ax2.set_xlim([-2, 2])
 ax2.set_xlabel('Eruption time in mins')
 ax2.set_ylabel('Waiting time to next eruption')
 ax2.set_title('Visualization of clustered data', y=1.02)
 ax2.set_aspect('equal')
 plt.tight_layout()
 plt.suptitle(f'Silhouette analysis using k = {k}',
fontsize=16, fontweight='semibold', y=1.05);

如上面的图所示，n_clusters=2的平均silhouette 得分最高，约为0.75左右，所有的聚类都高于均值，这说明它实际上是一个很好的选择。此外，silhouette 图的厚度可以显示出每个聚类有多大。从图中可以看出，第一组的样本数量几乎是第二组的两倍。然而，当我们将n_clusters增加到3和4时，平均silhouette 分数分别显著下降到0.48和0.39左右。此外，silhouette 图的厚度开始出现较大的波动。底线是:好的n_cluster具有远高于0.5 silhouette 平均得分，并且所有的clusters都高于平均值。

缺点

如果clusters具有球形形状，则Kmeans算法可以很好地捕获数据结构。它总是试图在质心周围构造一个漂亮的球形。这意味着，在群集具有复杂几何形状的那一刻，kmeans在聚类数据方面做得很差。我们将举例说明kmeans表现不佳的三种情况。

首先，kmeans算法不允许远离彼此的数据点共享同一个cluster，即使它们显然属于同一个cluster。下面是两条不同水平线上的数据点示例，说明了kmeans如何尝试将每条水平线的一半数据点组合在一起。Python代码如下：

# Create horizantal data
X = np.tile(np.linspace(1, 5, 20), 2)
y = np.repeat(np.array([2, 4]), 20)
df = np.c_[X, y]
 
km = KMeans(n_clusters=2)
km.fit(df)
labels = km.predict(df)
centroids = km.cluster_centers_
 
fig, ax = plt.subplots(figsize=(6, 6))
plt.scatter(X, y, c=labels)
plt.xlim([0, 6])
plt.ylim([0, 6])
plt.text(5.1, 4, 'A', color='red')
plt.text(5.1, 2, 'B', color='red')
plt.text(2.8, 4.1, 'C', color='red')
ax.set_aspect('equal')

Kmeans认为点'B'更接近点'A'而不是点'C'，因为它们具有非球形形状。因此，点'A'和'B'将位于同一cluster中，但点'C'将位于不同的cluster中。请注意，Single Linkage层次聚类方法是正确的，因为它不会分离类似的点）。

其次，我们将从具有不同均值和标准差的多元正态分布生成数据。因此，我们将有3组数据，其中每组由不同的多元正态分布（不同的平均值/标准偏差）生成。一组将拥有比其他两组更多的数据点。接下来，我们将对K = 3的数据运行kmeans，看看它是否能够正确地聚类数据。为了使比较更容易，我将首先根据它来自的分布绘制彩色数据。然后我将绘制相同的数据，但现在根据已分配的群集进行着色。Python代码如下：

# Create data from three different multivariate distributions
X_1 = np.random.multivariate_normal(mean=[4, 0], cov=[[1, 0], [0, 1]], size=75)
X_2 = np.random.multivariate_normal(mean=[6, 6], cov=[[2, 0], [0, 2]], size=250)
X_3 = np.random.multivariate_normal(mean=[1, 5], cov=[[1, 0], [0, 2]], size=20)
df = np.concatenate([X_1, X_2, X_3])
 
# Run kmeans
km = KMeans(n_clusters=3)
km.fit(df)
labels = km.predict(df)
centroids = km.cluster_centers_
 
# Plot the data
fig, ax = plt.subplots(1, 2, figsize=(10, 10))
ax[0].scatter(X_1[:, 0], X_1[:, 1])
ax[0].scatter(X_2[:, 0], X_2[:, 1])
ax[0].scatter(X_3[:, 0], X_3[:, 1])
ax[0].set_aspect('equal')
ax[1].scatter(df[:, 0], df[:, 1], c=labels)
ax[1].scatter(centroids[:, 0], centroids[:, 1], marker='o',
 c="white", alpha=1, s=200, edgecolor='k')
for i, c in enumerate(centroids):
 ax[1].scatter(c[0], c[1], marker='$%d$' % i, s=50, alpha=1, edgecolor='r')
ax[1].set_aspect('equal')
plt.tight_layout()

看来kmeans不能正确地计算出聚类。因为它试图最小化聚类内部的变化，所以它给大clustesr比小clusters更大的权重。换句话说，较小clusters中的数据点可能会远离质心，以便更多地关注较大的clusters。

最后，我们将生成具有复杂几何形状的数据，例如彼此之间的卫星和圆圈，并在这两个数据集中测试kmeans。Python代码如下：

# Cricles
X1 = make_circles(factor=0.5, noise=0.05, n_samples=1500)
 
# Moons
X2 = make_moons(n_samples=1500, noise=0.05)
 
fig, ax = plt.subplots(1, 2)
for i, X in enumerate([X1, X2]):
 fig.set_size_inches(18, 7)
 km = KMeans(n_clusters=2)
 km.fit(X[0])
 labels = km.predict(X[0])
 centroids = km.cluster_centers_
 
 ax[i].scatter(X[0][:, 0], X[0][:, 1], c=labels)
 ax[i].scatter(centroids[0, 0], centroids[0, 1], marker='*', s=400, c='r')
 ax[i].scatter(centroids[1, 0], centroids[1, 1], marker='+', s=300, c='green')
plt.suptitle('Simulated data', y=1.05, fontsize=22, fontweight='semibold')
plt.tight_layout()

正如预期的那样，kmeans无法找出两个数据集的正确聚类。但是，如果我们使用kernel 方法，我们可以帮助kmeans完美地聚类这些数据集。我们的想法是转换为更高维度的表示，使数据线性可分（与我们在SVM中使用的想法相同）。不同类型的算法在这种情况下非常有效，例如SpectralClustering，见下文：

# Cricles
X1 = make_circles(factor=0.5, noise=0.05, n_samples=1500)
 
# Moons
X2 = make_moons(n_samples=1500, noise=0.05)
 
fig, ax = plt.subplots(1, 2)
for i, X in enumerate([X1, X2]):
 fig.set_size_inches(18, 7)
 sp = SpectralClustering(n_clusters=2, affinity='nearest_neighbors')
 sp.fit(X[0])
 labels = sp.labels_
 ax[i].scatter(X[0][:, 0], X[0][:, 1], c=labels)
plt.suptitle('Simulated data', y=1.05, fontsize=22, fontweight='semibold')
plt.tight_layout();

结论

Kmeans聚类是最流行的机器学习聚类算法之一，通常是从业者在解决聚类任务时应用的第一件事，以了解数据集的结构。kmeans的目标是将数据点分成不同的、不重叠的子组。当clusters有一种球形形状时，它做得很好。然而，当clusters的几何形状偏离球面形状时，它就会受到影响。