KMeans Algorithm

K-Means 알고리즘은 주어진 데이터를 미리 정해진 개수 \(K\)의 클러스터로 분할하는 대표적인 클러스터링 알고리즘이다. 각 클러스터는 하나의 중심점(centroid)으로 대표되며, 알고리즘은 각 데이터를 가장 가까운 중심점에 할당하는 방식으로 데이터를 분할한다.

먼저 \(K\)개의 초기 클러스터 중심을 설정한다. 각 데이터는 가장 가까운 중심을 갖는 클러스터에 할당되고, 각 클러스터의 중심은 그 클러스터에 속한 데이터의 평균(무게중심)으로 다시 계산된다. 이후 새로운 중심을 기준으로 데이터를 다시 분할하고, 클러스터 중심을 다시 계산하는 과정을 반복한다. 이러한 반복은 클러스터의 중심 또는 데이터의 소속이 더 이상 변하지 않을 때까지 계속된다.

알고리즘의 수행 과정은 다음과 같다.

1. \(K\)개의 초기 클러스터 중심을 설정한다.
2. 각 데이터를 가장 가까운 클러스터에 할당한다.
3. 각 클러스터에 속한 데이터의 평균을 계산하여 새로운 클러스터 중심으로 설정한다.
4. 클러스터 중심 또는 데이터의 소속이 더 이상 변하지 않을 때까지 2와 3을 반복한다.

K-Means는 각 데이터와 자신이 속한 클러스터 중심 사이의 거리 제곱의 합

\begin{align}J=\sum_{k=1}^{K}\sum_{j\in C_k}\left| \mathbf{x}_j - \boldsymbol{\mu}_k \right|^2\end{align}

을 최소화하는 방향으로 클러스터를 분할한다. 여기서 \(C_k\)는 \(k\)번째 클러스터에 속한 데이터의 집합이고, \(\boldsymbol{\mu}_k\)는 해당 클러스터의 중심이다. 각 반복 과정에서 이 비용 함수는 감소하거나 유지되므로 알고리즘은 항상 종료되지만, 일반적으로는 전역 최적해가 아닌 국소 최적해를 찾는다.

K-Means의 결과는 초기 클러스터 중심의 선택에 영향을 받으므로 서로 다른 초기값에서 서로 다른 분할 결과가 얻어질 수 있다. 또한 클러스터의 개수 \(K\)를 미리 지정해야 한다는 단점이 있다. 실제 데이터에서는 적절한 \(K\)를 알기 어려운 경우가 많으며, 잘못된 \(K\)의 선택은 클러스터링 성능을 크게 저하시킬 수 있다. 이러한 문제를 해결하기 위해 클러스터의 분할과 병합을 반복하면서 클러스터의 개수를 자동으로 조정하는 ISODATA 알고리즘이 제안되었다.

typedef double * vector_t ;
double kmeans_label(
             const vector_t *obs,             // observations;
             const int n_obs,                   /* in # of vectors */            
             const vector_t *mean,           // mean vectors;
             const int n_mean,                /* # of mean vectors */
             int* label,                           // cluster labeling;
             const int veclen)                 // dimension of feature;
{
    double sqerr=0;
    for (int i = 0; i < n_obs; i++) {
        vector_t c = obs[i];
        /* Get an estimate of best distance (min_dist) and (min_idx) */
        int min_idx = label[i];
        vector_t m = mean[min_idx];
        double min_sqdist = 0 ;
        for (int l = 0; l < veclen; l++) {
            double t = m[l] - c[l]; min_sqdist += t * t;
        }
        for (int j = 0; j < n_mean; j++) {
            vector_t m = mean[j];
            // save time by checking:: dist < min_sqdist=previous min_val;
            double sqdist = 0;
            for (l = 0; (l < veclen) && (sqdist < min_sqdist); l++) {
                double t = m[l] - c[l]; sqdist += t * t;
            }
            if (sqdist < min_sqdist) {
                min_sqdist = sqdist; min_idx = j;
            }
        }
        label[i] = min_idx;
        sqerr += min_sqdist;
    }
    return sqerr ;
};

int kmeans_update(
              const vector_t *obs, //observations;
              const int n_obs,       //# of observations;     
              vector_t *mean,       // mean vectors;
              const int n_mean,    // # of mean vectors;
              int *label,               // cluster labelling;
              const int veclen      // dimension of features;
              ) {
    std::vector<int> cnt(n_mean);
    for (int i = 0; i < n_mean; i++)
        for (int l = 0; l < veclen; l++) mean[i][l] = 0.0;
        
    for (i = 0; i < n_obs; i++) {
        ASSERT((0 <= label[i]) && (label[i] < n_mean));  
        vector_t m = mean[label[i]];
        cnt[label[i]]++;  
        vector_t c = obs[i];
        if (c == NULL)
            TRACE("No observations for %u, but expected up through %u", i, n_obs-1);
        for (int l = 0; l < veclen; l++) m[l] += c[l];
    }
    for (i = 0; i < n_mean; i++) {
        int j = cnt[i];
        if (j != 0) for (int l = 0; l < veclen; l++) mean[i][l] /= (double) j;
        else {
            TRACE("Empty cluster %u\n", i);
            return FALSE;
        }
    }
    return TRUE;
}

double kmeans(
       const vector_t *obs,     // observations;
       const int n_obs,          /* # of observations */
       vector_t *mean,           /* initial set of means */
       const int n_mean,        /* # of means (should be k_mean?) */
       int *label,                   /* cluster labeling*/
       const int veclen,         /* vector length of means and corpus */    
       const double min_conv_ratio,
       const int max_iter
    ) {
    double p_sqerr = DBL_MAX;
    int ret;
    double sqerr = kmeans_label( obs, n_obs,mean, n_mean, &label[0], veclen);
    double conv_ratio = (p_sqerr - sqerr) / p_sqerr;
    for (int i = 0; (i < max_iter) && (conv_ratio > min_conv_ratio); i++) {
        TRACE("kmeans iter [%u] %e ...\n", i, conv_ratio);
        ret = kmeans_update(obs, n_obs, mean, n_mean, &label[0], veclen );
        if (ret != TRUE)
            return (double)ret;
        p_sqerr = sqerr;
        sqerr = kmeans_label(  obs, n_obs,mean, n_mean, &label[0], veclen);
        conv_ratio = (p_sqerr - sqerr) / p_sqerr;
    }
    TRACE("kmeans n_iter %u sqerr %e conv_ratio %e\n", i, sqerr, conv_ratio);   
    return sqerr;
}