k-meansを改良したk-means++を説明する

# はじめに
　先日、代表的なクラスタリング手法であるk-meansについて解説記事を書きました。
 
https://qiita.com/sasshi_i/items/2af0a76b95176437f923
 
しかし、k-meansには問題点があり、今回はその問題点に対応したk-means++という手法について説明します。

# k-meansの問題点
k-meansは最初のステップでクラスタの中心(セントロイド)をランダムにデータ点から選択します。ランダムにデータ点を選んでいることから以下のような問題が発生します。

* 初期値によって最適なクラスタリングが行われない可能性がある
* 初期値によって結果の収束までに時間がかかる可能性がある

# k-meansの問題点に対応したk-means++
## k-means++とは？
k-meansの初期値選択に改良を行ったクラスタリング手法です。初期のk個のクラスタはなるべく離れている方が良いという考えに基づいてセントロイドを配置します。

## k-means++のアルゴリズム
k-means++は以下の手順でクラスタリングを行います。
1. **初期クラスタ中心の設置**: データセットの中からランダムに1点選んでクラスタの中心とする
2. **距離の計算**: データセット内の各データと最近傍のクラスタ中心までの距離$D(x)$を計算する
3. **重み付き確率分布の計算**: データセット内の各データにおいて、重み付き確率分布$\frac{D(x)^2}{\sum D(x)^2}$を算出
4. **新しいクラスタ中心の設置**: 重み付き確率分布を用いて、新しいクラスタ中心を設置します。重み付き確率分布が大きい(既存のクラスタ中心から遠い)点ほどクラスタ中心に選ばれやすくなります
5. **繰り返し**:k個のクラスタ中心ができるまで、2~4を繰り返します
6. **k-meansの実施**: 設置されたk個の初期クラスタ中心を用いてk-meansを実施する

## ライブラリを使わずにk-means++を実装する
ここからはscikit-learnなどのライブラリを使わずにk-means++を実装したいと思います。

まずはデータの生成。
```python
import numpy as np
import random
from matplotlib import pyplot as plt

np.random.seed(0)
points1 = np.random.normal(5, 4, (50, 2))
points2 = np.random.normal(20, 4, (50, 2))
points3 = np.random.normal(30, 4, (50, 2))
points = np.r_[points1, points2, points3]
np.random.shuffle(points)

for point in points:
  plt.scatter(point[0], point[1], marker="o", color="blue")
plt.show()
```
![](https://storage.googleapis.com/zenn-user-upload/ad656452efc5-20240115.png)

k-means++の実装
```python
def calc_euclidean_distances(centroids, data):
  norms = []
  for centroid in centroids:
    norms.append([np.linalg.norm(centroid - data_point) for data_point in data])
  return np.array(norms).T
  
def k_means_pp(n_cluster, data, iter=1000):
  centroids = []

  # 初期クラスタ中心の設置
  init_centroid_index = random.randint(0, data.shape[0])
  centroids.append(data[init_centroid_index])

  # 選択した点を入力データから除外
  _data = np.delete(data, init_centroid_index, axis=0)

  # セントロイドとの距離の計算
  while len(centroids) < n_cluster:
    dists = []

    # 各データ点と最近傍のセントロイドとの距離を算出
    for data_point in _data:
      dists_for_point = []
      for centroid in centroids:
        dist = np.linalg.norm(centroid - data_point)
        dists_for_point.append(dist)
      dists.append(min(dists_for_point))

    # 重み付き確率分布の計算
    prob = dists/np.sum(dists)

    # 新しいクラスタ中心の設置
    new_c_id = np.random.choice([i for i in range(len(dists))], p=prob)
    centroids.append(_data[new_c_id])

    _data = np.delete(_data, new_c_id, axis=0)

  # k-meansの実施
  for _ in range(iter):
  # セントロイドとデータ点との距離を計算
    distances = calc_euclidean_distances(centroids, data)

    # セントロイドとの距離が最小となるクラスターに割り当て
    clusters = np.argmin(distances, axis=1)

    # クラスタごとに平均を算出
    for i_cluster in range(n_cluster):
      # 該当のクラスタに割り当てられたデータを抽出
      data_indices = [i for i, x in enumerate(clusters) if x == i_cluster]
      data_in_cluster = [data[i] for i in data_indices]

      # 平均値を算出し、セントロイドを更新
      centroids[i_cluster] = np.array(data_in_cluster).mean(axis=0)
  return clusters
```

データがうまくクラスタリングできるか確認。今回はクラスタ数3で検証します。
```python
clusters = k_means_pp(3, points)
for cluster, point in zip(clusters,points):
  if cluster == 0:
    plt.scatter(point[0], point[1], marker="o", color="blue")
  elif cluster == 1:
    plt.scatter(point[0], point[1], marker="o", color="green")
  else:
    plt.scatter(point[0], point[1], marker="o", color="red")
plt.show()
```
![](https://storage.googleapis.com/zenn-user-upload/d1de1892c36b-20240115.png)

# 終わりに
 k-meansでクラスタリングした際は、初期値によってうまくクラスタに分けられない時がありました。k-means++では、初期値によって上手くクラスタに分けられないという事象が明らかに減ったような気がします。
また、ライブラリを使わずに実装したことによってk-means++の仕組みを理解することができました。本記事が皆さんのお役に立つと嬉しいです。何か間違っているところがあれば、コメントなどでご指摘ください。

# 参考
https://github.com/scikit-learn/scikit-learn/blob/main/sklearn/cluster/_kmeans.py#L178

# 関連記事
https://qiita.com/sasshi_i/items/2af0a76b95176437f923