k-means clustering with custom distance function and k-means++

k-means는 Euclidean 공간에 주어진 $N$개의 데이터를 $k$개의 클러스터로 나누는 대표적인 비지도 클러스터링(군집화) 알고리즘으로 $k$ 개 클러스터 분산(sum of squared error)의 합이 최소가 되는 것을 목표로 한다. 본 포스트에서는 이 k-means 알고리즘의 수학적인 설명과 초기값 선택 문제를 다룰 것이다. 그리고 클러스터링 할 때 Euclidean 거리를 다른 특정한 거리로 적용하는 문제를 다루고자 한다. 먼저 k-means에 대해 살펴보자.

k-means

정의

데이터 $X = (\mathbf x_1, …, \mathbf x_N)$, $\mathbf x_i \in \mathbb R^d$에 대해서 $X$의 클러스터(또는 파티션)를 $S=${ $S_1, … S_k$} 라고 하자. 즉,

\[\bigcup_{i=1}^k S_i = X, \quad S_i \cap S_j = \emptyset, \forall i \neq j.\]

k-means는 클러스터 $S$의 분산 또는 within-cluster sum of squared(wcss)를 최소화하는 클러스터 $S$를 찾는 문제이다. 즉, 다음 목적함수를 최소화한다.

\[\underset{S}{\operatorname{argmin}} \sum_{i=1}^k \sum_{\mathbf x \in S_i} \| \mathbf x - \mu_i \|^2 \label{obj}\tag{1}\]

여기서 $\mu_i$는 클러스터 $S_i$의 샘플평균을 의미한다.

\[\mu_i:= \frac{1}{|S_i|} \sum_{\mathbf x \in S_i} \mathbf x\]

Remark

  • k-means는 NP-hard 문제로 글로벌 해를 찾는 것은 어렵기 때문에 heuristic algorithm으로 로컬 해를 찾는다.
  • 아이디어는 다음과 같다. $S=${ $S_1, … S_k$}가 (\ref{obj}) 목적함수의 글로벌 해라고 가정해보자. 그렇다면 모든 $i=1,…,k$에 대해서 $S_i$에 있는 샘플을 기준으로 Euclidean 거리가 가장 가까운 샘플평균은 $\mu_i$이다. 즉, $\mathbf x\in S_i$이면 다음 부등식을 만족한다.
\[\| \mathbf x- \mu_i \| \leq \| \mathbf x - \mu_j \| \, \forall j, 1\leq j \leq k\]
  • (\ref{obj}) 목적함수가 최소가 되는 클러스터인 경우 샘플들은 각 클러스터의 샘플평균과 가까운 점들의 모임으로 해석할 수 있다(물론 특정 샘플이 두 샘플평균과 거리가 동일한 거리에 있다면 유일하게 결정되지 않는다). 그래서 heuristic algorithm은 클러스터 $S=${ $S_1, … S_k$}를 샘플평균으로 매개화하여 글로벌 해 샘플평균 $\mu_1, …, \mu_k$을 찾는다.

Standard algorithm

  • 센트로이드(centroid)로 클러스터를 매개화(센트로이드에 가까이 있는 점들의 모임)하고 클러스터 내에 있는 샘플들과 센트로이드의 Euclidean 거리 제곱합이 최소가 되는 방식으로 센트로이드를 업데이트한다.
  • $k$개 초기 센트로이드 $m_1^{(1)},…, m_k^{(1)}$, $m_i^{(1)}\in \mathbb R^d$는 $N$개의 데이터 중에서 임의로(uniformly) $k$개로 선택한다(아래 k-means++는 초기 센트로이드를 최대한 멀리 떨어져있게 선택하는 방법이다). 그리고 아래 두 스텝을 반복한다.

  • (Assignment step) $t\geq 1$ 단계의 센트로이드와 Euclidean 거리 제곱이 가장 가까운 점들을 모아서 $S^{(t)}=${ $S_1^{(t)}, …, S_k^{(t)}$ }을 다음과 같이 정의한다.
\[S_i^{(t)}:= \{ \mathbf x \in X: \| \mathbf x- m_i^{(t)} \|^2 \leq \| \mathbf x - m_j^{(t)} \|^2 \, \forall j, 1\leq j \leq k \}.\]
  • (Update step) 위 단계에서 얻은 클러스터 $S^{(t)}=${ $S_1^{(t)}, … S_k^{(t)}$ }로부터 센트로이드를 샘플평균으로 업데이트한다.
\[m_i^{(t+1)}:= \frac{1}{|S_i^{(t)}|} \sum_{\mathbf x \in S_i^{(t)}} \mathbf x \quad (\text{sample mean of }S_i^{(t)})\]
Initial centroid & Assignment Update step
Assignment step Update step

k-means++

  • 일반적으로 k-means 알고리즘은 초기 센트로이드 선택이 클러스터링 결과에 큰 영향을 미친다. 초기 센트로이드가 모여 있는 경우 좋은 결과를 얻기 어렵다. 그래서 초기 센트로이드 $c_1,…, c_k$를 최대한 멀리 떨어져있는 점을 선택하기 위한 방법으로 Arthur, Vassilvitskii는 k-means++ 방법을 제시하였다.
    • 데이터 $X$에서 임의로 한 점을 선택해서 $c_1$으로 둔다.
    • 모든 샘플 $\mathbf x_i$와 $c_1$과의 거리를 계산해서 다음과 같은 확률로 다음 두 번째 센트로이드 $c_2$를 결정한다(쉽게 보면 $c_1$과 가장 먼 점이 선택될 확률이 높다).
    \[\frac{\| \mathbf x_i - c_1\|^2}{\sum_{j=1}^N \| \mathbf x_i - c_1\|^2}\]
    • $j\geq 3$에 대해서 다음을 수행하여 $c_j$ 을 결정한다. 각 샘플에서 각 센트로이드 까지 계산하고 가장 가까운 센트로이드에 할당한다. $i= 1, …, N$ 및 $p=1,…,j-1$에 대해서 다음 확률로 $c_j$를 선택한다. 이 방법으로 $j=k$가 될 때까지 $c_j$를 선택한다.
    \[\frac{\|\mathbf x_i - c_p\|^2}{\sum_{\mathbf x_j\in C_p} \|\mathbf x_j - c_p\|^2}\]
    • 여기서 $C_p$는 센트로이드 $c_p$에 가까운 점들을 모아놓은 집합이고 $\mathbf x_i$는 $C_p$에 속한다.
    • 예를 들어서, $j=3$ 인 경우 $c_1$과 가까이 있는 점들 $C_1$은 $c_1$과의 거리에 비례해서 선택될 확률이 결정되고 $c_2$와 가까이 있는 점들 $C_2$는 $c_1$과의 거리에 비례해서 선택될 확률이 결정된다. 즉, 기존에 선택된 $c_1$, $c_2$ 두 점과의 거리를 고려해서 그 다음 센트로이드 $c_3$를 결정한다. 이미 선택된 가장 가까운 센트로이드와 거리에 비례하는 확률로 후속 센트로이드가 선택된다고 보면 된다.
  • Python Scikit learn의 k-means는 k-means++로 초기 센트로이드를 선택하는 것이 기본(디폴트) 설정이다.

Python code

  • 아래 코드는 Scikit learn의 k-means 코드를 변형한 것이다.
import numpy as np
import matplotlib.pyplot as plt
from sklearn.cluster import KMeans
from sklearn import datasets
from mpl_toolkits.mplot3d import Axes3D


iris = datasets.load_iris()
X = iris.data
y = iris.target

# 디폴트 초기 센트로이드는 k-means++
# init = "random" 으로 설정 변경 시 임의의 선택으로 변경됨

estimator = KMeans(n_clusters=3, init="random")
estimator.fit(X)


fig = plt.figure(figsize=(6, 6))
ax = Axes3D(fig, rect=[0, 0, .95, 1], elev=48, azim=134)
labels = estimator.labels_

ax.scatter(X[:, 3], X[:, 0], X[:, 2],
           c=labels.astype(np.float), edgecolor='k')

ax.w_xaxis.set_ticklabels([])
ax.w_yaxis.set_ticklabels([])
ax.w_zaxis.set_ticklabels([])
ax.set_xlabel('Petal width')
ax.set_ylabel('Sepal length')
ax.set_zlabel('Petal length')
ax.set_title('k=3')
ax.dist = 12

fig = plt.figure(figsize=(6, 6))
ax = Axes3D(fig, rect=[0, 0, .95, 1], elev=48, azim=134)

for name, label in [('Setosa', 0),
                    ('Versicolour', 1),
                    ('Virginica', 2)]:
    ax.text3D(X[y == label, 3].mean(),
              X[y == label, 0].mean(),
              X[y == label, 2].mean() + 2, name,
              horizontalalignment='center',
              bbox=dict(alpha=.2, edgecolor='w', facecolor='w'))
# Reorder the labels to have colors matching the cluster results
y = np.choose(y, [1, 2, 0]).astype(np.float)
ax.scatter(X[:, 3], X[:, 0], X[:, 2], c=y, edgecolor='k')

ax.w_xaxis.set_ticklabels([])
ax.w_yaxis.set_ticklabels([])
ax.w_zaxis.set_ticklabels([])
ax.set_xlabel('Petal width')
ax.set_ylabel('Sepal length')
ax.set_zlabel('Petal length')
ax.set_title('Ground Truth')
ax.dist = 12

fig.show()
k-means(k=3) Ground truth

k-means 특정 거리로 일반화 하기

k-means의 알고리즘의 Euclidean 거리를 일반 거리인 $d(\cdot, \cdot)$로 일반화한다는 것이 무엇인지 살펴보자. 간단하게 k-means의 목적함수 (\ref{obj})의 Euclidean 거리만 $d(\cdot, \cdot)$으로 변경한 식인

\[\underset{S}{\operatorname{argmin}} \sum_{i=1}^k \sum_{\mathbf x \in S_i} d( \mathbf x, \mu_i)\]

을 최소화하는 문제로 생각할 수 있는데 이는 수학적으로 맞지 않다. 그 이유는 k-means의 목적함수 (\ref{obj})에 있는 샘플평균 $\mu_i$는 데이터 셋(클러스터) $S_i$에 있는 점들과 Euclidean 거리 제곱 합이 최소가 되는 점(센트로이드)이기 때문이다. 다시 말해서 데이터 셋 $S$ 의 샘플평균은 아래식($d$가 Euclidean 거리 제곱인 경우)을 만족하는 점이다.

\[\underset{\mathbf z \in \mathbb R^d}{\operatorname{argmin}} \sum_{\mathbf x\in S} \| \mathbf z - \mathbf x \|^2\]

데이터 셋 $S$와 거리 $d$에 대해서

\[\underset{\mathbf z \in \mathbb R^d}{\operatorname{argmin}} \sum_{\mathbf x\in S} d( \mathbf z, \mathbf x )\]

식을 만족하는 점을 거리 $d$에 대한 $S$의 센트로이드라고 하자. 그렇다면 거리 $d$에 대한 k-means는 다음 목적함수를 최소화 하는 것이다.

\[\underset{S}{\operatorname{argmin}} \sum_{i=1}^k \sum_{\mathbf x \in S_i} d( \mathbf x, \mathbf m_i)\]

여기서 $\mathbf m_i$는 거리 $d$에 대한 $S_i$의 센트로이드이다. 기본 k-means는 $d$가 Euclidean 거리의 제곱인 경우(샘플평균이 센트로이드)이고 $d$가 맨해튼거리이면 맨해튼거리에 대한 센트로이드는 중앙값(median)이다.

기본 k-means와 마찬가지로 $d$로 k-means로 일반화한 것도 NP-hard 문제이다. 그래서 k-means의 Standard Algorithm으로 클러스터 해를 찾는다. 즉, 다음과 같은 방법으로 해를 찾는다.

Standard algorithm of k-means implementation with custom distance

  • 센트로이드(centroid)로 클러스터를 매개화(센트로이드에 가까이 있는 점들의 모임)하고 클러스터 내에 있는 샘플들과 센트로이드의 일반 거리 $d$가 최소가 되는 방식으로 센트로이드를 업데이트한다.
  • $k$개 초기 센트로이드 $m_1^{(1)},…, m_k^{(1)}$, $m_i^{(1)}\in \mathbb R^d$는 $N$개의 데이터 중에서 임의로(uniformly) $k$개로 선택한다. 그리고 아래 두 스텝을 반복한다.

  • (Assignment step) $t\geq 1$ 단계의 센트로이드와 $d$ 거리가 가장 가까운 점들을 모아서 $S^{(t)}=${ $S_1^{(t)}, …, S_k^{(t)}$ }을 다음과 같이 정의한다.
\[S_i^{(t)}:= \{ \mathbf x \in X: d(\mathbf x, m_i^{(t)}) \leq d( \mathbf x, m_j^{(t)}) \, \forall j, 1\leq j \leq k \}.\]
  • (Update step) 위 단계에서 얻은 클러스터 $S^{(t)}=${ $S_1^{(t)}, … S_k^{(t)}$ }로부터 거리 $d$에 대한 센트로이드로 업데이트한다.
\[m_i^{(t+1)}:=\text{거리 } d \text{에 대한 }S_i\text{ 의 센트로이드}\]

Pyclustering for k-means implementation with custom distance

Pyclustering에서는 k-means를 일반 거리로 일반화하는 알고리즘을 제공한다. 그렇지만 일반적인 거리 $d$에 대해서 그 거리에 대한 센트로이드를 찾는 것은 어려운 일이다. 예를 들어서, $d$가 $L^p$, $p>1$ 거리일 때 이 거리에 대한 센트로이드를 찾는 것은 쉽지 않다.

Pyclustering에서 제공하는 일반화된 k-means 코드를 보자. (Assignment step) 에서는 센트로이드와 $d$ 거리가 가장 가까운 점들을 모으는 것은 맞는 방법이다. 하지만 (Update step) 에서는 거리 $d$에 대한 센트로이드로 업데이트 하는 것이 아니라 기본 k-means 처럼 샘플평균으로 업데이트한다.

Pyclustering의 업데이트 코드(거리와 관계없이)

def __update_centers(self):
    """!
    @brief Calculate centers of clusters in line with contained objects.
    
    @return (numpy.array) Updated centers.
    
    """
    
    dimension = self.__pointer_data.shape[1]
    centers = numpy.zeros((len(self.__clusters), dimension))
    
    for index in range(len(self.__clusters)):
        cluster_points = self.__pointer_data[self.__clusters[index], :]
        centers[index] = cluster_points.mean(axis=0)

    return numpy.array(centers)

위 방법이 수학적으로 makes sense 하지 않지만 Pyclustering은 k-means의 거리를 일반화하는 코드를 제공한다(그래도 일반화 해야 한다면). 그리고 잘 알려진 다양한 거리를 적용할 수 있다. 참조

IRIS 데이터에 대해 Cosine similarity로 k-means(Pyclustering)을 돌리는 코드는 다음과 같다.

# pip install pyclustering

import numpy as np
from numpy.linalg import norm
import matplotlib.pyplot as plt
from sklearn import datasets

from mpl_toolkits.mplot3d import Axes3D

from pyclustering.cluster.kmeans import kmeans
from pyclustering.cluster.center_initializer import kmeans_plusplus_initializer
from pyclustering.utils.metric import type_metric, distance_metric

iris = datasets.load_iris()
X = iris.data
y = iris.target

# 클러스터 개수
n_k = 3

# 일반화된 거리에 대한 k-means++
initial_centers = kmeans_plusplus_initializer(X, n_k).initialize()

# 거리가 Cosine similarity인 경우
def user_function(x, y):
    return np.dot(x,y) / (norm(x) * norm(y))

metric = distance_metric(type_metric.USER_DEFINED, func=user_function)

kmeans_instance = kmeans(X, initial_centers, metric = metric)
# run cluster analysis and obtain results
kmeans_instance.process()
clusters = kmeans_instance.get_clusters()

# 동일 클러스터에 있는 인덱스를 보여줌
print(clusters)
fig = plt.figure(figsize=(6, 6))
ax = Axes3D(fig, rect=[0, 0, .95, 1], elev=48, azim=134)
labels = kmeans_instance.predict(X)

ax.scatter(X[:, 3], X[:, 0], X[:, 2],
           c=labels.astype(np.float), edgecolor='k')

ax.w_xaxis.set_ticklabels([])
ax.w_yaxis.set_ticklabels([])
ax.w_zaxis.set_ticklabels([])
ax.set_xlabel('Petal width')
ax.set_ylabel('Sepal length')
ax.set_zlabel('Petal length')
ax.set_title('k=3')
ax.dist = 12
[[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49], [50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 145, 146, 147, 148, 149]]

Reference