Кластеризация

Все модели до сих пор, от ближайших соседей до глубоких сетей из предыдущей заметки, учились на парах $(x, y)$: вход и ответ, который он должен давать. Разметка — самая дорогая часть набора данных, и часто её нет вовсе: журнал пользовательских сессий, миллион неразмеченных изображений, векторы, которые выдаёт энкодер. У данных всё равно есть структура, и первый вопрос к ним — распадаются ли они на группы. Эта заметка об этом вопросе. Она показывает, почему «группе» нужно определение, прежде чем у вопроса появится ответ, разбирает самое употребительное определение, метод k-средних, как задачу оптимизации с простым алгоритмом, а затем переходит к определениям, которые k-средних выразить не может: плотные области, иерархии и графы.

Обозначения. $x_1, \dots, x_n \in \mathbb R^d$ — точки данных, $K$ — число кластеров, $c_i \in \{1, \dots, K\}$ — кластер точки $i$, $C_k = \{i : c_i = k\}$ — точки кластера $k$, а $n_k = \lvert C_k \rvert$ — их число. $\mu_1, \dots, \mu_K \in \mathbb R^d$ — центры кластеров (центроиды), $\lVert \cdot \rVert$ — евклидова норма.

Что такое кластер?#

Кластеризация (clustering) — это разбиение точек на группы, такое что точки внутри группы похожи, а точки из разных групп — нет. Это определение оставляет открытыми два выбора: что значит «похожи» и какой формы могут быть группы. Три распространённых ответа:

• Компактные группы. Каждый кластер — шар вокруг центра, и в хорошей кластеризации расстояния от точек до их центра малы. Это метод k-средних.
• Плотные области. Кластер — область, где точки расположены тесно, любой формы, отделённая от других кластеров разреженным пространством; изолированные точки не принадлежат ни одному кластеру. Это DBSCAN.
• Связная структура. Кластер — множество точек, связанных цепочками близких соседей, или часть графа сходства, из которой выходит мало рёбер. Это одиночная связь и спектральная кластеризация.

На хорошо разделённых круглых облаках все три определения согласны. На всём остальном они дают разные ответы, и ни один из них не ошибочен: они отвечают на разные вопросы. Клейнберг (2002) сделал это утверждение точным. Он искал функцию кластеризации от попарных расстояний с тремя свойствами: она не зависит от единицы измерения расстояний (масштабная инвариантность, scale invariance), любое разбиение получается для каких-то расстояний (полнота, richness), а уменьшение расстояний внутри кластеров при увеличении расстояний между ними не меняет ответа (согласованность, consistency). Его теорема 2.1 показывает, что ни одна функция не обладает всеми тремя свойствами, а теорема 2.2 — что одиночная связь с подходящим правилом остановки обеспечивает любые два из них. Кластеризации без определения не бывает.

Выбор расстояния важен не меньше выбора метода. Все методы ниже видят данные только через расстояния, поэтому единицы измерения признаков определяют ответ: умножение одного признака на 10 делает его главным в каждом расстоянии. На данных из четырёх облаков из раздела «Код» метод k-средних восстанавливает облака точно (скорректированный индекс Рэнда 1.0, определение ниже), а после умножения второго признака на 10 — лишь частично (0.736). Масштабирование и кодирование признаков, о которых заметка о данных и признаках, — часть кластеризации, а в высоких размерностях сами расстояния становятся менее информативными по причине, объяснённой в заметке о ближайших соседях и в разделе о сферическом слое многомерного нормального распределения.

k-средних: целевая функция#

Метод k-средних (k-means) выражает «компактность» числом. Для назначения $c$ и центров $\mu$ внутрикластерная сумма квадратов (within-cluster sum of squares) равна

$$J(c, \mu) = \sum_{i=1}^n \lVert x_i - \mu_{c_i} \rVert^2 = \sum_{k=1}^K \sum_{i \in C_k} \lVert x_i - \mu_k \rVert^2 ,$$

и метод k-средних ищет $c$ и $\mu$, которые её минимизируют. В scikit-learn $J$ называется инерцией (inertia). Два факта делают целевую функцию наглядной.

Лучший центр кластера — его среднее. При фиксированных точках $C_k$ функция $\mu \mapsto \sum_{i \in C_k} \lVert x_i - \mu \rVert^2$ — строго выпуклая квадратичная; её градиент $-2\sum_{i \in C_k} (x_i - \mu)$ обращается в ноль в среднем $\bar x_k = \frac{1}{n_k}\sum_{i \in C_k} x_i$. Разложение вокруг среднего показывает, сколько теряет любой другой центр:

$$\sum_{i \in C_k} \lVert x_i - \mu \rVert^2 = \sum_{i \in C_k} \lVert x_i - \bar x_k \rVert^2 + n_k \lVert \bar x_k - \mu \rVert^2 ,$$

потому что перекрёстный член $2(\bar x_k - \mu)^\top \sum_{i \in C_k}(x_i - \bar x_k)$ равен нулю. Это утверждение заметки об условном математическом ожидании о том, что среднее — лучшее постоянное предсказание при квадратичной ошибке, применённое к каждому кластеру.

Малый разброс внутри означает большой разброс между. Суммирование того же тождества по кластерам при $\mu = \bar x$, среднем всех данных, даёт

$$\underbrace{\sum_{i=1}^n \lVert x_i - \bar x \rVert^2}_{\text{полный}} = \underbrace{\sum_{k=1}^K \sum_{i \in C_k} \lVert x_i - \bar x_k \rVert^2}_{\text{внутри}} + \underbrace{\sum_{k=1}^K n_k \lVert \bar x_k - \bar x \rVert^2}_{\text{между}} ,$$

выборочный вариант закона полной дисперсии. Полный разброс задан данными, поэтому минимизировать разброс внутри кластеров — то же самое, что максимизировать разброс средних кластеров. Внутрикластерный член можно записать и без центров, $\sum_{i \in C_k} \lVert x_i - \bar x_k \rVert^2 = \frac{1}{2 n_k} \sum_{i, j \in C_k} \lVert x_i - x_j \rVert^2$: метод k-средних штрафует средний квадрат расстояния между точками одного кластера.

Точная минимизация $J$ NP-трудна, даже для двух кластеров, если размерность входит в размер входа (Aloise et al., 2009). На практике используют локальный поиск.

Алгоритм Ллойда — покоординатный спуск#

$J$ зависит от двух блоков переменных: дискретного назначения $c$ и непрерывных центров $\mu$. Каждый блок легко оптимизировать при фиксированном другом, и алгоритм Ллойда (Lloyd’s algorithm; Lloyd, 1982; название k-means восходит к MacQueen, 1967) чередует две точные минимизации:

1. Шаг назначения. При фиксированных $\mu$ функция $J$ — сумма отдельных слагаемых, по одному на точку, поэтому каждая точка независимо выбирает ближайший центр: $c_i \leftarrow \arg\min_k \lVert x_i - \mu_k \rVert^2$ (при равенстве точка остаётся в текущем кластере).
2. Шаг обновления. При фиксированном $c$ функция $J$ — сумма отдельных слагаемых, по одному на кластер, и каждое минимизируется средним кластера: $\mu_k \leftarrow \bar x_k$. (Кластер, потерявший все точки, сохраняет старый центр или перезапускается в далёкой точке.)

Каждый шаг минимизирует $J$ по одному блоку, поэтому ни один не может её увеличить:

$$J(c^{t+1}, \mu^{t+1}) \le J(c^{t+1}, \mu^{t}) \le J(c^{t}, \mu^{t}) .$$

Это покоординатный спуск (coordinate descent) по двум координатам. Он ещё и завершается. Назначений не больше $K^n$; когда шаг назначения меняет $c$, он строго уменьшает $J$ (точка переходит только к строго более близкому центру), поэтому ни одно назначение не повторяется, и алгоритм останавливается за конечное число шагов в точке, которую не может улучшить ни один из шагов. Эта точка — локальный минимум в этом покоординатном смысле, а не глобальный. Каждая итерация стоит $O(nKd)$ на вычисление расстояний, и на практике итераций нужно немного.

Разобранный пример. Пять точек на прямой, $0, 2, 3, 8, 10$, при $K = 2$ и начальных центрах $\mu = (0, 3)$.

итерация	центры	кластеры	$J$
0	$0,\ 3$	$\{0\}$, $\{2, 3, 8, 10\}$	$0 + 1 + 0 + 25 + 49 = 75$
1	$0,\ 5.75$	$\{0, 2\}$, $\{3, 8, 10\}$	$34.69$
2	$1,\ 7$	$\{0, 2, 3\}$, $\{8, 10\}$	$16$
3	$1.67,\ 9$	$\{0, 2, 3\}$, $\{8, 10\}$	$6.67$

На итерации 3 назначение больше не меняется, средние остаются на месте, и алгоритм останавливается с $J = 6.67$. В двух измерениях то же чередование выглядит так:

Рисунок

Алгоритм Ллойда: назначение и пересчёт центров

Шаг назначения разрезает плоскость на ячейки Вороного (Voronoi cells): точки, которые ближе к $\mu_k$, чем к $\mu_j$, образуют полуплоскость по одну сторону от серединного перпендикуляра к отрезку между $\mu_k$ и $\mu_j$, а ячейка — пересечение таких полуплоскостей, выпуклый многоугольник. Поэтому любая кластеризация методом k-средних — разбиение на выпуклые части с прямыми границами, какой бы формы ни были данные. Одно это предложение объясняет большинство неудач, описанных ниже.

Инициализация: k-means++#

Алгоритм Ллойда сходится к тому локальному минимуму, к которому ведёт спуск из начальной точки. На рисунке два из трёх начальных центров оказываются в одном облаке, и алгоритм справляется только потому, что один из них может дойти до пустого облака. При других стартах этого не происходит: на данных из четырёх облаков из раздела «Код» один запуск из четырёх случайных точек выборки достигает лучшего значения $J = 1127.4$ в 56% из 200 запусков, а среднее по запускам — $J = 4064.4$, потому что неудачные запуски сливают два облака и делят другое.

k-means++ (Arthur and Vassilvitskii, 2007) выбирает начальные центры так, чтобы они были разнесены. Первый центр — равномерно выбранная случайная точка выборки. Каждый следующий центр — точка выборки, выбранная с вероятностью, пропорциональной $D(x)^2$, квадрату расстояния от $x$ до ближайшего из уже выбранных центров. Вероятнее всего выбираются точки, далёкие от всех текущих центров, — скорее всего, они лежат в кластере без центра. Теорема 3.1 их статьи оценивает результат одной только инициализации, до каких-либо итераций Ллойда:

$$\mathbb E[J] \le 8(\ln K + 2)\, J_{\text{opt}} .$$

При $K = 4$ множитель около 27 — оценка грубая, но это гарантия, которой нет у равномерного старта; на практике инициализация оказывается гораздо ближе к оптимуму. На тех же данных каждый из 200 одиночных стартов k-means++ достигает $J = 1127.4$. (k-means++ в scikit-learn берёт для каждого нового центра несколько кандидатов и оставляет тот, что сильнее всего уменьшает $J$, — жадный вариант.) Обычная практика добавляет перезапуски: алгоритм запускают с нескольких начальных значений генератора и оставляют наименьшее $J$.

Разобранный пример. Для точек $0, 2, 3, 8, 10$ и первого центра в $0$ квадраты расстояний равны $0, 4, 9, 64, 100$, их сумма $177$. Второй центр — $10$ с вероятностью $100/177 = 0.565$ и $8$ с вероятностью $0.362$: с вероятностью $0.93$ он попадает в дальнюю группу, тогда как равномерный выбор попал бы туда с вероятностью $0.4$.

Сколько кластеров?#

$J$ сама по себе не может выбрать $K$. Оптимальное значение $J$ не возрастает с ростом $K$, поскольку ещё один центр всегда можно поставить на уже существующий, и обращается в 0 при $K = n$, когда у каждой точки свой кластер. Две распространённые эвристики определяют $K$ по данным.

«Локоть» (elbow). Строим лучшее значение $J$ в зависимости от $K$ и ищем точку, после которой добавление кластеров перестаёт окупаться. На данных из четырёх облаков $J = 18363,\ 5241,\ 1127,\ 1020,\ 921$ для $K = 2, \dots, 6$: падение от 3 к 4 велико, падения после 4 малы, и локоть приходится на $K = 4$. На реальных данных кривая часто изгибается плавно, и положение локтя — дело вкуса.

Силуэт (silhouette; Rousseeuw, 1987) оценивает каждую точку по тому, насколько она ближе к своему кластеру, чем к следующему. Пусть $a(i)$ — среднее расстояние от $x_i$ до остальных точек её кластера, а $b(i)$ — наименьшее среднее расстояние от $x_i$ до точек другого кластера. Тогда

$$s(i) = \frac{b(i) - a(i)}{\max\big(a(i), b(i)\big)} \in [-1, 1] ,$$

причём $s(i) = 0$ для точки, которая одна в своём кластере. Значения около 1 означают, что точка глубоко внутри своего кластера, около 0 — на границе, отрицательные — что она, вероятно, попала не в тот кластер. Среднее по всем точкам оценивает кластеризацию, и его можно сравнивать при разных $K$. Для пяти точек $0, 1, 2, 6, 7$ с кластерами $\{0, 1, 2\}$ и $\{6, 7\}$ у точки $2$ получаем $a = (2 + 1)/2 = 1.5$, $b = (4 + 5)/2 = 4.5$ и $s = (4.5 - 1.5)/4.5 = 0.67$. На данных из четырёх облаков средний силуэт равен $0.588, 0.762, 0.811, 0.680, 0.577$ для $K = 2, \dots, 6$, наибольший при $K = 4$. Силуэт использует все попарные расстояния, $O(n^2)$, и, как и $J$, предпочитает выпуклые, хорошо разделённые кластеры: он наследует представление о кластере из метода k-средних. Вероятностная альтернатива, штрафованное правдоподобие, появится вместе со смесями распределений в следующей заметке.

Где k-средних ошибается#

Целевая функция предпочитает кластеры круглые, с похожим разбросом и похожего размера, а геометрия Вороного ничего другого выразить не может.

• Невыпуклые кластеры. Две переплетённые полулуны нельзя разделить прямой, поэтому k-средних разрезает каждую пополам. На 400 точках двух лун из раздела «Код» скорректированный индекс Рэнда с настоящими лунами равен 0.274.
• Разный разброс. Граница между двумя кластерами — серединный перпендикуляр к отрезку между их центрами, посередине между ними, каким бы ни был разброс. Если плотный кластер (стандартное отклонение 0.5) находится в четырёх единицах от широкого (стандартное отклонение 2), внешние точки широкого кластера, лежащие ближе к центру плотного, назначаются плотному; скорректированный индекс Рэнда равен 0.561. Гауссовская смесь, в которой у каждого кластера своя ковариационная матрица, достигает на тех же данных 0.902; ей посвящена следующая заметка.
• Разные размеры и вытянутые кластеры. Разбиение большого кластера может уменьшить $J$ сильнее, чем сохранение маленького, и вытянутый кластер по той же причине разрезается поперёк длинной оси.
• Выбросы. Среднее смещается к далёким точкам, и одна далёкая точка может забрать себе целый центр. Замена квадратов расстояний модулями, а средних — медианами (k-медиан), или центров — точками выборки (k-медоидов) делает метод устойчивее ценой больших вычислений.

Рисунок

Две луны: k-средних против DBSCAN

DBSCAN: кластеры как плотные области#

DBSCAN (Ester et al., 1996) определяет кластеры через плотность и требует двух параметров: радиуса $\varepsilon$ и числа $\text{minPts}$. Обозначим через $N_\varepsilon(x)$ точки выборки на расстоянии не больше $\varepsilon$ от $x$, включая саму $x$.

• Точка — основная (core point), если $\lvert N_\varepsilon(x)\rvert \ge \text{minPts}$: она находится в плотной области.
• Две основные точки на расстоянии не больше $\varepsilon$ друг от друга принадлежат одному кластеру, а по цепочке — и все основные точки, соединённые путём из таких шагов. Кластер — максимальное множество основных точек, связанных таким образом, вместе с неосновными точками на расстоянии не больше $\varepsilon$ от какой-либо из них, его граничными точками (border points).
• Все остальные точки — шум (noise) и не принадлежат ни одному кластеру.

Алгоритм наращивает кластеры поиском в ширину из непосещённых основных точек. Кластеры могут быть любой формы, потому что цепочка $\varepsilon$-шагов может следовать кривой; число кластеров — выход, а не вход; выбросы помечаются как выбросы, а не втискиваются в какой-то кластер. Граничная точка, лежащая в пределах $\varepsilon$ от двух кластеров, достаётся тому, который дотянется до неё первым, — это единственная зависимость от порядка. Каждый шаг требует запроса окрестности: при вычислении всех попарных расстояний стоимость $O(n^2)$, а пространственный индекс делает её гораздо меньше на данных малой размерности.

Разобранный пример. Возьмём на прямой точки $0, 0.5, 1, 1.5, 5, 9, 9.4, 9.8$ с $\varepsilon = 0.6$ и $\text{minPts} = 3$. Окрестность точки $0.5$ — $\{0, 0.5, 1\}$, окрестность точки $1$ — $\{0.5, 1, 1.5\}$, поэтому обе основные; у $0$ и $1.5$ в пределах $0.6$ только по две точки, и это граничные точки того же кластера. Аналогично $9.4$ — основная точка с граничными точками $9$ и $9.8$. У точки $5$ нет соседей: это шум. Два кластера и один выброс, причём $K$ не задавали.

Выбор параметров. В исходной статье для двумерных данных берут $\text{minPts} = 4$, а $\varepsilon$ определяют по отсортированному графику $k$-расстояний (sorted k-distance graph): расстоянию от каждой точки до её $k$-го ближайшего соседа, $k = \text{minPts}$, отсортированному по убыванию. У точек шума $k$-расстояния большие, у точек кластеров — маленькие, и $\varepsilon$ ставят в изломе между ними (раздел 4.2). Метод чувствителен к $\varepsilon$: на лунах из раздела «Код» $\varepsilon = 0.1$ даёт четыре кластера и 10 точек шума, $\varepsilon = 0.2$ точно восстанавливает две луны, а $\varepsilon = 0.3$ сливает всё в один кластер.

Одна плотность для всех кластеров. Единственное $\varepsilon$ задаёт один порог плотности. Когда один кластер гораздо плотнее другого, радиус, достаточно малый, чтобы отделить плотный кластер от соседей, объявляет разреженный кластер шумом. HDBSCAN (Campello, Moulavi and Sander, 2013) рассматривает все радиусы сразу: кластеризации при убывающем $\varepsilon$ образуют иерархию, и он оставляет самые устойчивые кластеры — те, что сохраняются в широком диапазоне радиусов.

Иерархическая кластеризация и способы связи#

Дополнительный материал: целая иерархия кластеризаций и способы связи, которые её строят. Остальная часть заметки на него не опирается.

Часто единственно верного $K$ нет: документы образуют темы, темы делятся на подтемы. Агломеративная кластеризация (agglomerative clustering) строит всю иерархию. Она начинает с $n$ кластеров из одной точки и раз за разом сливает два ближайших кластера, пока не останется один. $n - 1$ слияний образуют дерево, дендрограмму (dendrogram), высоты слияний в которой показывают, насколько далеки были сливаемые кластеры; разрез дерева на некоторой высоте или на уровне, где остаётся $K$ кластеров, даёт плоскую кластеризацию.

«Ближайшим кластерам» нужно расстояние между множествами точек — способ связи (linkage):

Связь	Расстояние между кластерами $A$ и $B$	Поведение
одиночная (single)	$\min_{a \in A, b \in B} \lVert a - b \rVert$	идёт по цепочкам; находит изогнутые кластеры, но один мостик из шума сливает два кластера
полная (complete)	$\max_{a \in A, b \in B} \lVert a - b \rVert$	компактные кластеры похожего диаметра
средняя (average)	среднее $\lVert a - b \rVert$ по всем парам	между одиночной и полной
Уорда (Ward)	прирост $J$ из-за слияния	целевая функция k-средних, жадно

Связь Уорда (Ward, 1963) сливает ту пару, которая меньше всего увеличивает внутрикластерную сумму квадратов. По разложению из раздела о k-средних слияние $A$ и $B$ увеличивает её на

$$\Delta(A, B) = \frac{n_A n_B}{n_A + n_B}\, \lVert \bar x_A - \bar x_B \rVert^2 .$$

Для $A = \{0, 2, 3\}$ и $B = \{8, 10\}$ из примера с k-средних $\Delta = \frac{3 \cdot 2}{5}(9 - 1.67)^2 = 64.53$ — ровно разность между разбросом всех пяти точек вокруг их среднего, $71.2$, и внутрикластерным разбросом $6.67$.

Одиночная связь — замаскированное минимальное остовное дерево: её кластеризация на $K$ групп получается удалением $K - 1$ самых длинных рёбер минимального остовного дерева точек (Gower and Ross, 1969). Поэтому она идёт по цепочкам близких соседей вдоль лун (скорректированный индекс Рэнда 1.0 в разделе «Код»), тогда как полная, средняя связь и связь Уорда разрезают их, как k-средних (0.285, 0.329 и 0.239). На круглых облаках все четыре восстанавливают кластеры. Цена иерархии — матрица попарных расстояний, $O(n^2)$ памяти, и от $O(n^2)$ до $O(n^3)$ времени в зависимости от способа связи, что ограничивает её десятками тысяч точек.

Спектральная кластеризация#

Дополнительный материал: кластеры как слабо связанные части графа сходства, найденные по собственным векторам его лапласиана; предполагается знание собственных значений и векторов симметричной матрицы. Остальная часть заметки на него не опирается.

Одиночная связь смотрит на цепочки; спектральная кластеризация (spectral clustering) смотрит на весь граф близких соседей и находит части, между которыми мало рёбер. Построим на точках граф сходства с весами $w_{ij} \ge 0$: например, $w_{ij} = 1$, если $x_j$ входит в число $m$ ближайших соседей $x_i$ или наоборот, или $w_{ij} = \exp(-\lVert x_i - x_j \rVert^2 / 2\sigma^2)$. Со степенями $d_i = \sum_j w_{ij}$ и $D = \operatorname{diag}(d_1, \dots, d_n)$ лапласиан графа (graph Laplacian) равен $L = D - W$. Для любого вектора $f \in \mathbb R^n$, по одному значению на точку,

$$f^\top L f = \frac12 \sum_{i,j=1}^n w_{ij} (f_i - f_j)^2 ,$$

и эта величина мала, когда $f$ почти постоянен на сильно связанных точках. Значит, $L$ симметричен и положительно полуопределён, для постоянного вектора $Lf = 0$, а вектор $f$ с $f^\top L f = 0$ обязан быть постоянным вдоль каждого ребра. Отсюда ключевой факт (von Luxburg, 2007, предложение 2): кратность собственного значения 0 матрицы $L$ равна числу компонент связности графа, а индикаторные векторы компонент порождают соответствующее собственное подпространство.

Если граф распадается ровно на $K$ компонент, то $K$ собственных векторов с собственным значением 0 — это, с точностью до поворота, индикаторы компонент; в координатах, которые они дают каждой точке, все точки одной компоненты совпадают. Если компоненты соединены несколькими слабыми рёбрами, собственные значения не нулевые, а малые, и собственные векторы почти кусочно-постоянны. Именно это и использует алгоритм:

1. Построить граф сходства и его лапласиан.
2. Вычислить собственные векторы $u_1, \dots, u_K$ для $K$ наименьших собственных значений.
3. Дать точке $i$ новые координаты $(u_{1i}, \dots, u_{Ki}) \in \mathbb R^K$ и запустить на них k-средних.

В новых координатах луны превращаются в два плотных облака, и k-средних их разделяет. На лунах из раздела «Код» граф 10 ближайших соседей имеет ровно две компоненты, у лапласиана два нулевых собственных значения, за которыми следует $0.02$, и k-средних по двум собственным векторам восстанавливает луны (скорректированный индекс Рэнда 1.0). Метод можно вывести и как релаксацию разреза графа: найти $K$ групп с малым весом рёбер между ними относительно их размера — NP-трудная задача (RatioCut, нормализованный разрез), а замена дискретных индикаторных векторов вещественными превращает её в задачу на собственные значения выше; нормализованным вариантам $L$ соответствуют нормализованный разрез Shi and Malik (2000) и алгоритм Ng, Jordan and Weiss (2001) (von Luxburg, раздел 5). Стоимость — спектральное разложение, $O(n^3)$ для плотной матрицы и гораздо меньше для разреженного графа ближайших соседей, а ответ зависит от того, как построен граф.

Оценка кластеризации#

Без разметки кластеризацию можно судить только по тому определению, из которого она построена. $J$ сравнивает решения с одинаковым $K$; силуэт сравнивает разные $K$, но наследует представление о кластере из k-средних: на лунах из раздела «Код» он даёт настоящему разбиению $0.33$, а разбиению k-средних, которое режет обе луны, $0.49$. Практическая проверка — устойчивость (stability): кластеризовать несколько бутстрэп-выборок или случайных подмножеств и измерить, насколько результаты согласуются. Устойчивость — свидетельство, а не доказательство. Структура, которая меняется при каждой новой выборке, подозрительна, но и настоящие кластеры, которые перекрываются или выражены слабо, бывают неустойчивыми, а устойчивый ответ может навязать сам метод: при $K = 2$ k-средних почти одинаково на каждой новой выборке режет одно вытянутое гауссовское облако поперёк длинной оси на две половины. Полезна ли кластеризация, зависит ещё от задачи и от представления, в котором считались расстояния.

При наличии разметки части данных, скрытой от кластеризации (сами точки кластеризуются, скрыты только их метки), кластеризацию можно сравнить с ней, но не по доле правильных ответов: номера кластеров произвольны, и кластеризация $\{1, 1, 2, 2\}$ совпадает с $\{2, 2, 1, 1\}$. Меры на основе подсчёта пар обходятся без сопоставления номеров. По всем $\binom n2$ парам точек считаем пары, на которых два разбиения согласны (оба помещают пару вместе или оба — порознь). Доля согласий — индекс Рэнда (Rand index). У него есть недостаток: два случайных разбиения согласны на многих парах случайно, в основном «порознь», поэтому его значения скапливаются около 1. Скорректированный индекс Рэнда (adjusted Rand index, ARI) Hubert and Arabie (1985) вычитает согласие, ожидаемое случайно. Пусть $n_{ij}$ — число точек истинного класса $i$ в кластере $j$, а $a_i$ и $b_j$ — суммы по строкам и столбцам. Тогда

$$\text{ARI} = \frac{\sum_{ij} \binom{n_{ij}}{2} - E}{\frac12\big[\sum_i \binom{a_i}{2} + \sum_j \binom{b_j}{2}\big] - E}, \qquad E = \frac{\sum_i \binom{a_i}{2} \sum_j \binom{b_j}{2}}{\binom n2} .$$

$E$ — ожидаемое число пар, находящихся вместе в обоих разбиениях, если перемешать метки кластеров при сохранении их размеров: гипергеометрическая модель случайности. ARI равен 1 для одинаковых разбиений, около 0 для несвязанных и может быть отрицательным. Точкам, которые DBSCAN называет шумом, нужно соглашение. Раздел «Код» передаёт в ARI метку шума $-1$ так, будто это ещё один кластер: при $\varepsilon = 0.1$ это даёт $0.731$, а оценка только по 390 кластеризованным точкам даёт $0.769$ — большее число, которое не учитывает 10 отброшенных точек. В сообщаемом ARI нужно указывать, какой из вариантов посчитан.

Разобранный пример. Шесть точек с истинными классами $AAABBB$ и кластерами $112222$. Таблица $n_{ij}$ имеет строки $(2, 1)$ и $(0, 3)$, поэтому $\sum \binom{n_{ij}}{2} = 1 + 0 + 0 + 3 = 4$; классы дают $\binom32 + \binom32 = 6$, кластеры — $\binom22 + \binom42 = 7$, и $\binom62 = 15$. Тогда $E = 6 \cdot 7 / 15 = 2.8$ и $\text{ARI} = (4 - 2.8)/(6.5 - 2.8) = 0.32$. Индекс Рэнда для той же пары разбиений равен $10/15 = 0.67$ — звучит гораздо лучше, чем есть на самом деле.

Разметка — лишь одна структура из многих. Изображения цифр можно группировать по цифре, по автору или по толщине штриха; низкий ARI относительно меток цифр может означать другую и всё же полезную группировку. Внешний индекс измеряет согласие с выбранной структурой, а не качество вообще.

Где k-средних появляется снова#

Метод k-средних — стандартный способ превратить непрерывные векторы в небольшой набор дискретных символов, и курс генеративных моделей несколько раз использует его именно так. VQ-VAE заменяет каждый вектор энкодера ближайшей записью выученной кодовой книги, и её обновление скользящим средним — это k-средних, вычисляемый онлайн, с той же проблемой: центр, который не выбирает ни одна точка, никогда не сдвигается. Нейросетевые аудиокодеки складывают такие кодовые книги в остаточное векторное квантование, где каждая ступень кластеризует то, что оставили предыдущие. ImageGPT превращает пиксели в токены с помощью палитры из 512 цветов, найденной методом k-средних, а кривая точности–полноты для генеративных моделей сравнивает гистограммы по кластерам k-средних для эмбеддингов настоящих и сгенерированных данных. В векторном поиске среди приближённых индексов, упомянутых в заметке о семантическом поиске, есть произведённое квантование (product quantization; Jégou, Douze and Schmid, 2011), которое делит каждый вектор на части и заменяет каждую часть ближайшим из набора центроидов k-средних.

Что изменилось#

	Обучение с учителем	k-средних	DBSCAN	Иерархическая	Спектральная
Вход	пары $(x, y)$	точки, $K$	точки, $\varepsilon$, minPts	точки, способ связи	граф сходства, $K$
Кластер — это		ячейка Вороного вокруг среднего	плотная связная область	поддерево дендрограммы	слабо связанная часть графа
Формы		выпуклые	любые	любые (одиночная) или компактные	любые
Выбросы		тянут центры	помечаются как шум	одноточечные ветви (одиночная)	зависит от графа
Стоимость		$O(nKd)$ на итерацию	$O(n^2)$ без индекса	$O(n^2)$ памяти	собственные векторы матрицы $n \times n$
Критерий оценки	ошибка на тесте	$J$, силуэт	ARI по разметке, устойчивость	дендрограмма, ARI	ARI по разметке, устойчивость

Ограничения#

• Ответ зависит от определения. Каждый метод закладывает своё представление о кластере, и ни один метод не прав на всех данных; выбор — часть моделирования.
• Расстояния должны быть осмысленными. Масштабирование и кодирование признаков меняют любой результат, а в высоких размерностях расстояния концентрируются; кластеризация эмбеддинга, выученного для этой цели, или нескольких главных компонент часто работает лучше кластеризации исходных признаков.
• Гиперпараметры без целевой переменной. $K$, $\varepsilon$, minPts и способ построения графа нельзя подобрать по ошибке на валидации, потому что её нет; эвристики для их выбора закладывают те же предположения, что и сами методы.
• Жёсткие назначения. Каждый метод выше помещает точку ровно в один кластер (или ни в один), без меры уверенности. Точка на границе двух кластеров обрабатывается так же, как точка в центре.
• Локальные оптимумы и стоимость. k-средних находит локальный минимум NP-трудной задачи; иерархическая и спектральная кластеризация масштабируются квадратично или хуже по $n$.

Что нужно запомнить#

• Кластеризации нужны определение кластера и расстояние; разные определения дают разные, одинаково допустимые ответы (теорема Клейнберга о невозможности).
• Метод k-средних минимизирует $J = \sum_i \lVert x_i - \mu_{c_i} \rVert^2$. Алгоритм Ллойда чередует назначение ближайшему центру и пересчёт средних; это покоординатный спуск, $J$ никогда не возрастает, и алгоритм останавливается в локальном минимуме.
• k-means++ выбирает начальные центры с вероятностью, пропорциональной $D(x)^2$, и гарантирует $\mathbb E[J] \le 8(\ln K + 2) J_{\text{opt}}$; его стоит использовать с несколькими перезапусками.
• $K$ выбирают по «локтю» кривой $J$ или по среднему силуэту $s = (b - a)/\max(a, b)$; одна $J$ всегда предпочитает больше кластеров.
• Ячейки k-средних выпуклы и делят пространство посередине между центрами, поэтому метод ошибается на изогнутых кластерах, при разном разбросе и при выбросах.
• DBSCAN: у основных точек не меньше minPts соседей в пределах $\varepsilon$; кластеры — цепочки основных точек вместе с их граничными точками; остальное — шум.
• Агломеративная кластеризация строит дендрограмму; одиночная связь идёт по цепочкам (минимальное остовное дерево), связь Уорда жадно минимизирует $J$.
• Спектральная кластеризация вкладывает точки с помощью нижних собственных векторов $L = D - W$, нулевые собственные значения которой считают компоненты связности, и запускает там k-средних.
• С разметкой кластеризацию оценивают скорректированным индексом Рэнда, который поправляет согласие пар на случайность; без разметки — силуэтом и устойчивостью, помня о том, что они предполагают.

Код#

Алгоритм Ллойда, инициализация, выбор $K$, случаи неудач и другие методы на лунах, со scikit-learn:

import numpy as np
from scipy.sparse.csgraph import laplacian
from sklearn.datasets import make_blobs, make_moons
from sklearn.cluster import KMeans, DBSCAN, AgglomerativeClustering
from sklearn.mixture import GaussianMixture
from sklearn.neighbors import kneighbors_graph
from sklearn.metrics import adjusted_rand_score as ari, silhouette_score

X, y = make_blobs(n_samples=600, centers=4, cluster_std=1.0, random_state=7)

# алгоритм Ллойда вручную из четырёх случайных точек выборки: J никогда не возрастает
rng = np.random.default_rng(0)
mu = X[rng.choice(len(X), 4, replace=False)]
for t in range(4):
    d = ((X[:, None, :] - mu[None]) ** 2).sum(axis=2)          # квадраты расстояний, n x K
    c = d.argmin(axis=1)                                         # шаг назначения
    print(t, round(d[np.arange(len(X)), c].sum(), 1))
    mu = np.array([X[c == k].mean(axis=0) for k in range(4)])   # шаг обновления: средние кластеров
# 0 22942.3 | 1 2406.5 | 2 1127.4 | 3 1127.4

# один случайный старт против одного старта k-means++, по 200 seed: среднее J и доля запусков с лучшим J
for init in ('random', 'k-means++'):
    J = np.array([KMeans(4, init=init, n_init=1, random_state=s).fit(X).inertia_ for s in range(200)])
    print(init, round(J.mean(), 1), (J < 1.001 * J.min()).mean())
# random 4064.4 0.56 | k-means++ 1127.4 1.0

# локоть и силуэт для K = 2..7
for K in range(2, 8):
    km = KMeans(K, n_init=10, random_state=0).fit(X)
    print(K, round(km.inertia_, 1), round(silhouette_score(X, km.labels_), 3))
# 2 18363.1 0.588 | 3 5241.2 0.762 | 4 1127.4 0.811 | 5 1019.9 0.68 | 6 921.1 0.577 | 7 813.9 0.446

# ответ зависит от единиц измерения: второй признак умножен на 10
km = lambda Z, K: KMeans(K, n_init=10, random_state=0).fit_predict(Z)
print(ari(y, km(X, 4)), round(ari(y, km(X * [1, 10], 4)), 3))       # 1.0 0.736

# разный разброс: k-средних против гауссовской смеси (следующая заметка)
Xu, yu = make_blobs(n_samples=[200, 200], centers=[[0, 0], [4, 0]], cluster_std=[0.5, 2.0], random_state=2)
print(round(ari(yu, km(Xu, 2)), 3), round(ari(yu, GaussianMixture(2, random_state=0).fit(Xu).predict(Xu)), 3))   # 0.561 0.902

# две луны: k-средних, DBSCAN, четыре способа связи
Xm, ym = make_moons(n_samples=400, noise=0.06, random_state=0)
print(round(ari(ym, km(Xm, 2)), 3))                                  # 0.274
print(round(silhouette_score(Xm, ym), 2), round(silhouette_score(Xm, km(Xm, 2)), 2))   # 0.33 0.49: силуэт предпочитает k-средних
for eps in (0.1, 0.2, 0.3):
    lab = DBSCAN(eps=eps, min_samples=5).fit_predict(Xm)
    print(eps, round(ari(ym, lab), 3), lab.max() + 1, (lab == -1).sum())   # ARI, кластеры, точки шума
# 0.1 0.731 4 10 | 0.2 1.0 2 0 | 0.3 0.0 1 0
lab = DBSCAN(eps=0.1, min_samples=5).fit_predict(Xm); kept = lab != -1
print(round(ari(ym[kept], lab[kept]), 3), kept.sum())               # 0.769 390: без точек шума
for link in ('single', 'complete', 'average', 'ward'):
    agg = lambda Z, K: AgglomerativeClustering(K, linkage=link).fit_predict(Z)
    print(link, round(ari(ym, agg(Xm, 2)), 3), round(ari(y, agg(X, 4)), 3))
# single 1.0 1.0 | complete 0.285 1.0 | average 0.329 1.0 | ward 0.239 1.0

# спектральная кластеризация вручную: граф 10 ближайших соседей, его лапласиан, k-средних по двум собственным векторам
W = kneighbors_graph(Xm, 10).toarray()
W = np.maximum(W, W.T)                                               # ребро, если одна из точек — сосед другой
lam, V = np.linalg.eigh(laplacian(W))                                # L = D - W
print(lam[:3].round(4), round(ari(ym, km(V[:, :2], 2)), 3))          # [0. 0. 0.02] 1.0: две компоненты, два нуля

Дальше#

Все методы этой заметки дают жёсткие ответы: точка принадлежит одному кластеру, и ничто не говорит, насколько это надёжно и насколько широким может быть каждый кластер. В частности, k-средних ведёт себя так, будто все кластеры круглые и одинаково широкие. Следующая заметка делает кластеризацию вероятностной. Каждый кластер становится гауссианой со своим средним и своей ковариационной матрицей, каждая точка получает вероятность принадлежать каждому кластеру, а параметры подбираются методом максимального правдоподобия с помощью EM-алгоритма, предельным случаем которого окажется алгоритм Ллойда.

Литература#

• J. Kleinberg. An Impossibility Theorem for Clustering. NeurIPS, 2002. Раздел 2 (три свойства, теорема 2.1, теорема 2.2).
• S. P. Lloyd. Least squares quantization in PCM. IEEE Transactions on Information Theory, 28(2), 129–137, 1982.
• J. MacQueen. Some methods for classification and analysis of multivariate observations. Proceedings of the Fifth Berkeley Symposium on Mathematical Statistics and Probability, 1, 281–297, 1967.
• D. Aloise, A. Deshpande, P. Hansen, P. Popat. NP-hardness of Euclidean sum-of-squares clustering. Machine Learning, 75(2), 245–248, 2009.
• D. Arthur, S. Vassilvitskii. k-means++: The Advantages of Careful Seeding. SODA, 2007. Раздел 2.2 (инициализация), теорема 3.1.
• P. J. Rousseeuw. Silhouettes: a graphical aid to the interpretation and validation of cluster analysis. Journal of Computational and Applied Mathematics, 20, 53–65, 1987.
• M. Ester, H.-P. Kriegel, J. Sander, X. Xu. A density-based algorithm for discovering clusters in large spatial databases with noise. KDD, 226–231, 1996. Раздел 3 (определения), раздел 4.2 (выбор Eps и MinPts).
• R. J. G. B. Campello, D. Moulavi, J. Sander. Density-Based Clustering Based on Hierarchical Density Estimates. PAKDD, LNCS 7819, 160–172, 2013.
• J. H. Ward. Hierarchical Grouping to Optimize an Objective Function. Journal of the American Statistical Association, 58(301), 236–244, 1963.
• J. C. Gower, G. J. S. Ross. Minimum Spanning Trees and Single Linkage Cluster Analysis. Applied Statistics, 18(1), 54–64, 1969.
• U. von Luxburg. A Tutorial on Spectral Clustering. Statistics and Computing, 17(4), 2007. Раздел 3 (предложения 1 и 2), раздел 4 (алгоритмы), раздел 5 (взгляд через разрезы графа).
• J. Shi, J. Malik. Normalized cuts and image segmentation. IEEE TPAMI, 22(8), 888–905, 2000.
• A. Y. Ng, M. I. Jordan, Y. Weiss. On Spectral Clustering: Analysis and an algorithm. NeurIPS, 2001.
• L. Hubert, P. Arabie. Comparing partitions. Journal of Classification, 2(1), 193–218, 1985.
• H. Jégou, M. Douze, C. Schmid. Product Quantization for Nearest Neighbor Search. IEEE TPAMI, 33(1), 117–128, 2011.
• T. Hastie, R. Tibshirani, J. Friedman. The Elements of Statistical Learning, 2nd ed. Springer, 2009. Раздел 14.3 (кластерный анализ).
• C. M. Bishop. Pattern Recognition and Machine Learning. Springer, 2006. Раздел 9.1 (кластеризация методом k-средних).
• E. Sokolov et al. Машинное обучение (HSE), конспекты лекций: лекция 12, обучение без учителя (2021); лекция 17, кластеризация (2025).