Участник:Ян-Мартин Тамм/Строгий алгоритм С средних

Строгий алгоритм С средних (K-means, Hard C-Means)
Последовательный алгоритм
Последовательная сложность	[math]\displaystyle{ O(nkdi) }[/math]
Объём входных данных	[math]\displaystyle{ nd+1 }[/math]
Объём выходных данных	[math]\displaystyle{ nk }[/math]
Параллельный алгоритм
Высота ярусно-параллельной формы	[math]\displaystyle{ O(log_2(k)log_2(d)i) }[/math]
Ширина ярусно-параллельной формы	[math]\displaystyle{ O(nkd) }[/math]

Описание алгоритма

Строгий алгоритм C средних, более известный как K-средних, решает задачу разделения векторов в многомерном пространстве на заданное количество кластеров. При этом каждый вектор принадлежит только одному кластеру, поэтому алгоритм называется строгим, в то время как нечеткий алгоритм C средних (Fuzzy C-means, Soft K-means) определяет степень принадлежности вектора каждому из кластеров.
Данная задача является NP-трудной, тем не менее описываемый алгоритм сходится к локальному оптимуму. Глобальный оптимум не гарантируется и зависит от первоначального распределения центров кластеров и их количества.
Общая идея алгоритма такова: фиксируем число кластеров [math]\displaystyle{ k }[/math], произвольным образом выбираем центры кластеров, после чего повторяем следующие два действия, пока не будет выполнено условие остановки:

1. Распределяем каждый вектор в ближайший кластер
2. Пересчитываем значения центров как среднее значение координат векторов в каждом кластере

Минимизируется сумма расстояний векторов от центров их кластеров. Остановка производится если либо значение данной функции становится ниже некоторого порога, либо мало изменяется по сравнению с предыдущей итерацией.

Математическое описание алгоритма

Дано множество d-мерных векторов  [math]\displaystyle{ U=\{u_1, u_2, ..., u_n\} }[/math] и число кластеров [math]\displaystyle{ k \in \mathbb{N} }[/math]. Требуется разбить множество [math]\displaystyle{ U }[/math] на [math]\displaystyle{ k }[/math] непересекающихся непустых кластеров [math]\displaystyle{ C_1, C_2, ..., C_k }[/math] таких, что [math]\displaystyle{ {\bigcup \limits _{i = 1}^k C_i} = U }[/math].
Для этого минимизируется функция [math]\displaystyle{ J = \sum\limits_{i=1}^k J_i = \sum\limits_{i=1}^k \sum\limits_{ u \in C_i} \|u-c_i \|^2 }[/math], где [math]\displaystyle{ c_i }[/math]— центр масс векторов кластера [math]\displaystyle{ C_i }[/math].
Произвольным образом инициализируем центры кластеров [math]\displaystyle{ c_1, c_2, ..., c_k }[/math], после чего начинается основная часть алгоритма.

• Каждый вектор [math]\displaystyle{ u_l }[/math] определяем в ближайший к нему кластер [math]\displaystyle{ C_i }[/math], то есть такой, что [math]\displaystyle{ \|u_l-c_i\|^2 \leqslant \|u_l-c_j\|^2 \;\; \forall i \neq j; \;\; (i,j=1,2,\ldots,k, \; l = 1,2,\ldots,n) }[/math].
Хранить информацию о принадлежности векторов кластерам будем в матрице [math]\displaystyle{ M }[/math] следующего вида:

[math]\displaystyle{ m_{ij} = \begin{cases} 1, & u_j \in C_i \\ 0, & \text{otherwise} \end{cases} }[/math]

• Считаем значение функции [math]\displaystyle{ J }[/math] и прекращаем работу алгоритма, если это значение меньше заданного порога [math]\displaystyle{ \varepsilon }[/math], либо значение мало изменилось после текущей итерации.
• Обновляем значения центров кластеров

  [math]\displaystyle{ c_i = \frac{1}{|C_i|} \sum\limits_{u \in C_i} u \;\; (i=1,2,\ldots,k) }[/math]

• Возвращаемся к распределению векторов по кластерам.

Вычислительное ядро алгоритма

Вычислительным ядром алгоритма является весь алгоритм, а именно два основных его шага:

• Распределение векторов по кластерам
• Пересчёт центров кластеров

Макроструктура алгоритма

Как уже говорилось раньше, первоначальное распределение векторов по кластерам может быть произвольным, например, можно случайным образом выбрать [math]\displaystyle{ k }[/math] векторов из множества всех векторов [math]\displaystyle{ U }[/math].
Далее начинается цикл, его первая операция — распределение векторов по кластерам. В ней для каждого вектора считается расстояние от него, до каждого кластера [math]\displaystyle{ \rho = \big \| u - c \big \|^2=\sum_{i=1}^d (u_i-c_i)^2 }[/math], после чего из них выбирается наименьшее и заносится в матрицу.
Подсчёт функции [math]\displaystyle{ J }[/math] можно проводить одновременно с первым шагом, так как в ней суммируются те же самые расстояния.
Пересчёт центров кластеров — сумма всех векторов кластера, делённая на количество объектов в кластере. [math]\displaystyle{ c = \frac{1}{|C|} \sum\limits_{u \in C} u }[/math]

Схема реализации последовательного алгоритма

def kmeans(points, k, cutoff):
    
    # Pick out k random points to use as our initial centroids
    initial = random.sample(points, k)
    
    # Create k clusters using those centroids
    clusters = [Cluster([p]) for p in initial]
    
    # Loop through the dataset until the clusters stabilize
    while True:
        # Create a list of lists to hold the points in each cluster
        lists = [ [] for c in clusters]

        # For every point in the dataset ...
        for p in points:
            # Get the distance between that point and the centroid of the first cluster.
            smallest_distance = getDistance(p, clusters[0].centroid)
        
            # Set the cluster this point belongs to
            clusterIndex = 0
        
            # For the remainder of the clusters ...
            for i in range(clusterCount - 1):
                # calculate the distance of that point to each other cluster's centroid
                distance = getDistance(p, clusters[i+1].centroid)
                # If it's closer to that cluster's centroid update what we
                # think the smallest distance is, and set the point to belong to that cluster
                if distance < smallest_distance:
                    smallest_distance = distance
                    clusterIndex = i+1
            lists[clusterIndex].append(p)
        
        # Set our total_shift to zero for this iteration
        total_shift = 0.0
        
        # As many times as there are clusters ...
        for i in range(k):
            # Calculate how far the centroid moved in this iteration
            shift = clusters[i].update(lists[i])
            # Keep track of the total moves from all cluster centroid updates
            total_shift += shift
        
        # If the centroids have stopped moving much, we're done!
        if total_shift < cutoff:
            break
    return clusters

Последовательная сложность алгоритма

Распределение по кластерам — для каждого из [math]\displaystyle{ n }[/math] векторов и каждого из [math]\displaystyle{ k }[/math] кластеров посчитать расстояние, сложив [math]\displaystyle{ d }[/math] раз значения, на вычисление каждого из которых нужно 2 операции — [math]\displaystyle{ O(nkd) }[/math]. Пересчёт центров — [math]\displaystyle{ n }[/math] раз сложить и [math]\displaystyle{ k }[/math] раз поделить —  [math]\displaystyle{ O(n+k) }[/math]. На завершение алгоритма понадобится [math]\displaystyle{ i }[/math] итераций. Итого сложность [math]\displaystyle{ O(nkdi) }[/math].

Ресурс параллелизма алгоритма

Самой трудоёмкой операцией является распределение по кластерам. Так, для рассчёта расстояний, можно провести суммирование [math]\displaystyle{ d }[/math] элементов за [math]\displaystyle{ log_2(d) }[/math] операций. Расстояние от данного вектора до всех кластеров можно найти за одну операцию, после чего нужно найти наименьшее значение из [math]\displaystyle{ k }[/math] элементов за [math]\displaystyle{ log_2(k) }[/math] операций. И всё это нужно сделать для [math]\displaystyle{ n }[/math] векторов, что можно сделать параллельно за одну операцию. Для сходимости всё так же потребуется [math]\displaystyle{ i }[/math] итераций, общая сложность получается [math]\displaystyle{ O(log_2(k)log_2(d)i) }[/math].

Входные и выходные данные алгоритма

Входные данные:

• [math]\displaystyle{ U=\{u_1, u_2, ..., u_n\} \subset \mathbb{R}^d }[/math] — множество векторов
• [math]\displaystyle{ k \in \mathbb{N} }[/math] — количество кластеров
• [math]\displaystyle{ \varepsilon }[/math] — порог остановки не обязательно.

Объём входных данных: [math]\displaystyle{ nd+1 }[/math].

Выходные данные:

• [math]\displaystyle{ M }[/math] — матрица принадлежности векторов кластерам размера [math]\displaystyle{ nk }[/math]. При этом данная матрица содержит только [math]\displaystyle{ n }[/math] единиц, остальные значения ноль. Возможен вывод в другом виде, например, список списков.

Объём выходных данных: [math]\displaystyle{ nk }[/math].

Свойства алгоритма

• Алгоритм является устойчивым — малая ошибка в начальных данных или при вычислениях очень незначительно повлияют на результат (скорее всего совсем не повлияют).