Кластеризация: k-means, DBSCAN, иерархическая, GMM

Все алгоритмы, которые мы разбирали до сих пор — от линейной регрессии до градиентного бустинга — учились на парах (x, y). Был правильный ответ, была функция потерь, измеряющая расхождение с ним, и был честный способ проверить себя: отложенная выборка.

Кластеризация устроена принципиально иначе. Есть только X — матрица n × d. Нужно разбить строки на группы так, чтобы «похожие оказались вместе». И тут начинается самое интересное: никто не сказал, что значит «похожие», и никто не скажет, правильно ли вы разбили.

Возьмём колоду карт и попросим трёх человек разложить её на группы. Первый разложит по мастям (4 группы), второй — по достоинству (13 групп), третий — на «картинки и числа» (2 группы). Все трое правы. Данные одни, а «естественная» структура зависит от того, какое понятие сходства вы неявно принесли с собой. Это не недостаток кластеризации, а её суть: алгоритм не находит структуру в данных, он находит структуру, соответствующую заложенной в него модели кластера. Ваша работа — выбрать модель, которая совпадает с бизнес-смыслом задачи.

Эта статья — про четыре разных ответа на вопрос «что такое кластер»:

Формальная постановка и почему «лучшей» кластеризации не бывает

Дано множество объектов X = {x₁, …, xₙ} и функция расстояния d(xᵢ, xⱼ). Требуется найти отображение C: X → {1, …, k} (жёсткая кластеризация) или C: X → Δᵏ — распределение по кластерам (мягкая кластеризация).

Отсутствие целевой переменной означает, что критерий качества приходится придумывать самому — и любой такой критерий будет одновременно определением кластера. Это не философия: Джон Клейнберг доказал теорему о невозможности кластеризации (An Impossibility Theorem for Clustering, NIPS 2002). Он выписал три естественных свойства, которых мы хотели бы от функции кластеризации:

1. Масштабная инвариантность — умножение всех расстояний на α > 0 не меняет результат.
2. Богатство — при подходящей метрике достижимо любое разбиение множества.
3. Согласованность — если сжать расстояния внутри кластеров и растянуть между кластерами, результат не изменится.

Ни одна функция кластеризации не удовлетворяет всем трём одновременно. Практический вывод: каждый алгоритм жертвует чем-то из списка. k-means с фиксированным k не обладает богатством. Single-linkage без порога отсечения не масштабно-инвариантен в части выбора уровня. Это освобождает от поиска «правильного алгоритма» и переводит вопрос в плоскость «какой компромисс подходит моей задаче».

Расстояние — половина решения

Прежде чем выбирать алгоритм, выберите метрику. Практически все алгоритмы этой статьи (кроме GMM, который работает с ковариациями) видят данные только через d(·,·).

• Евклидово ‖xᵢ − xⱼ‖₂ — по умолчанию; требует сопоставимых масштабов признаков.
• Манхэттенское ‖·‖₁ — устойчивее к выбросам в отдельных координатах.
• Косинусное 1 − cos(xᵢ, xⱼ) — стандарт для текстовых эмбеддингов и TF-IDF, где важно направление, а не длина вектора.
• Жаккар — для бинарных/множественных признаков (корзина покупок, набор тегов).
• Gower — смешанные типы (числовые + категориальные), редко реализовано, но иногда единственный вариант.

Первая по частоте ошибка новичка: подать k-means сырые колонки возраст (0–100) и доход (0–500000). Евклидово расстояние будет на 99.99% определяться доходом, возраст просто не существует для алгоритма. Масштабирование обязательно — детали в статье про данные и признаки.

Вторая ошибка — игнорировать проклятие размерности. В высокой размерности отношение (d_max − d_min) / d_min стремится к нулю: все точки становятся примерно одинаково далеки друг от друга, и понятие «ближайший сосед» теряет смысл (Beyer et al., When Is “Nearest Neighbor” Meaningful?, 1999). Поэтому на эмбеддингах размерности 768 сначала делают снижение размерности, а потом кластеризуют.

Карта алгоритмов

Дальше подробно разбираем четыре опорных семейства; остальные — вариации на их темы.

k-means: минимизация внутрикластерной дисперсии

Целевая функция

k-means ищет разбиение на k кластеров, минимизирующее сумму квадратов расстояний до центров (within-cluster sum of squares, WCSS, в sklearn — inertia_):

J(C, μ) = Σᵢ₌₁ⁿ ‖xᵢ − μ_{C(i)}‖²  →  min по C и μ

Ключевое наблюдение: при фиксированном разбиении оптимальный центр — это среднее точек кластера (отсюда название), а при фиксированных центрах оптимальное назначение — ближайший центр. Задача в целом NP-трудна даже для k = 2 в общей размерности, но каждый из двух подшагов решается тривиально. Это подсказывает алгоритм покоординатного спуска — алгоритм Ллойда (1957, опубликован в 1982).

Алгоритм Ллойда

Псевдокод:

KMEANS(X, k, max_iter, tol):
    μ ← KMEANSPP_INIT(X, k)
    повторять до max_iter:
        для каждой точки xᵢ:
            c[i] ← argmin_j ‖xᵢ − μⱼ‖²          # O(n·k·d)
        для каждого кластера j:
            μ'ⱼ ← среднее {xᵢ : c[i] = j}        # O(n·d)
            если кластер пуст: μ'ⱼ ← самая удалённая от своего центра точка
        если max_j ‖μ'ⱼ − μⱼ‖ < tol: стоп
        μ ← μ'
    вернуть c, μ

Почему сходится. Оба шага не увеличивают J: переназначение точки к ближайшему центру не может увеличить её вклад; пересчёт центра как среднего минимизирует сумму квадратов внутри кластера. Значит J монотонно не возрастает, а число возможных разбиений конечно — следовательно, алгоритм останавливается за конечное число итераций. Но останавливается в локальном минимуме, зависящем от инициализации, и этот минимум может быть сколь угодно хуже глобального.

k-means++ — почему инициализация решает

Случайный выбор k стартовых точек регулярно даёт вырожденные решения: два центра попадают в один плотный кластер, а два далёких кластера сливаются. k-means++ (Arthur & Vassilvitskii, 2007) выбирает центры последовательно, с вероятностью, пропорциональной квадрату расстояния до ближайшего уже выбранного центра (D²-sampling). Авторы доказали, что уже до запуска итераций Ллойда ожидаемое значение J не хуже 8(ln k + 2) от оптимума — логарифмическая гарантия аппроксимации взамен полного отсутствия гарантий.

В sklearn это поведение по умолчанию (init="k-means++", n_init="auto"). Отключать его без причины не нужно.

Реализация с нуля

import numpy as np


def kmeanspp_init(X: np.ndarray, k: int, rng: np.random.Generator) -> np.ndarray:
    """D^2-sampling: центры выбираются тем вероятнее, чем дальше они от уже выбранных."""
    n = X.shape[0]
    centers = np.empty((k, X.shape[1]), dtype=X.dtype)
    centers[0] = X[rng.integers(n)]
    # d2[i] — квадрат расстояния от x_i до ближайшего выбранного центра
    d2 = ((X - centers[0]) ** 2).sum(axis=1)
    for j in range(1, k):
        probs = d2 / d2.sum()
        centers[j] = X[rng.choice(n, p=probs)]
        d2 = np.minimum(d2, ((X - centers[j]) ** 2).sum(axis=1))
    return centers


def kmeans(X: np.ndarray, k: int, max_iter: int = 300, tol: float = 1e-4, seed: int = 0):
    """Алгоритм Ллойда. Возвращает (метки, центры, инерцию)."""
    rng = np.random.default_rng(seed)
    mu = kmeanspp_init(X, k, rng)
    labels = np.zeros(X.shape[0], dtype=np.int64)

    x_sq = (X ** 2).sum(axis=1, keepdims=True)

    for _ in range(max_iter):
        # ‖x − μ‖² = ‖x‖² − 2·x·μ + ‖μ‖² — одно матричное умножение вместо двойного цикла
        d2 = x_sq - 2 * (X @ mu.T) + (mu ** 2).sum(axis=1)   # (n, k)
        new_labels = d2.argmin(axis=1)
        own_d2 = d2[np.arange(X.shape[0]), new_labels]       # расстояние до своего центра

        new_mu = np.empty_like(mu)
        for j in range(k):
            mask = new_labels == j
            if mask.any():
                new_mu[j] = X[mask].mean(axis=0)
            else:
                # пустой кластер — переносим центр в худше всего описанную точку
                new_mu[j] = X[own_d2.argmax()]

        shift = np.linalg.norm(new_mu - mu, axis=1).max()
        mu, labels = new_mu, new_labels
        if shift < tol:
            break

    inertia = ((X - mu[labels]) ** 2).sum()
    return labels, mu, float(inertia)

Трюк с раскрытием квадрата нормы важен на практике: он превращает вычисление матрицы расстояний в одно матричное умножение, которое BLAS выполняет на порядок быстрее наивного двойного цикла.

Сложность и масштабирование

Для больших n есть MiniBatchKMeans: центры обновляются по случайным подвыборкам размера b, стоимость итерации падает до O(b · k · d). Качество (inertia) обычно на 1–3% хуже, скорость — в 10–100 раз выше (Sculley, Web-Scale K-Means Clustering, WWW 2010). Это рабочий выбор для десятков миллионов объектов.

Где k-means ломается

Слева — фундаментальное ограничение. k-means строит разбиение Вороного: граница между любыми двумя кластерами — гиперплоскость, перпендикулярная отрезку между их центрами. Значит все кластеры получаются выпуклыми. Кольцо внутри кольца, полумесяцы, вытянутые «колбасы» — всё это k-means разрежет геометрически, а не по смыслу.

Полный список слабых мест:

1. Только выпуклые, примерно сферические кластеры одинакового радиуса. Дисперсия штрафуется квадратично, поэтому большие и вытянутые кластеры алгоритм стремится разбить.
2. Чувствительность к масштабу признаков (см. выше) и к их корреляции: k-means не умеет учитывать ковариацию — это умеет GMM.
3. Чувствительность к выбросам. Один объект в 100 σ утаскивает центр на себя. Лечится либо чисткой данных, либо переходом на k-medoids (центром служит реальный объект, минимизирующий сумму расстояний), либо MiniBatchKMeans с усечёнными признаками.
4. k задаётся руками.
5. Кластеры примерно равного размера. Если реальные группы 95% / 5%, k-means охотно разрежет большую пополам, а маленькую растворит.

Как выбрать k

Единственно верного ответа нет — есть набор диагностик, которые вместе с доменным смыслом дают решение.

Метод локтя. Строим inertia(k). Функция монотонно убывает (при k = n она равна нулю), но в точке «истинного» k темп убывания резко замедляется. Проблема: излом часто размыт, а на реальных данных его может не быть вовсе.

Силуэт (Rousseeuw, 1987). Для точки i: a(i) — среднее расстояние до своих, b(i) — минимальное среднее расстояние до чужого кластера. Тогда

s(i) = (b(i) − a(i)) / max(a(i), b(i)) ∈ [−1, 1]

Значение около 1 — точка уверенно в своём кластере; около 0 — на границе; отрицательное — вероятно, попала не туда. Средний силуэт по выборке — интегральная оценка. Стоимость O(n²·d), для больших выборок считают на подвыборке (sample_size в sklearn).

Gap statistic (Tibshirani, Walther, Hastie, 2001) сравнивает log(inertia) на ваших данных с тем же на равномерном шуме в том же боксе. Выбирается наименьшее k, при котором разрыв достаточно велик. Единственный из перечисленных методов, умеющий сказать «кластеров нет вообще, k = 1».

Стабильность. Кластеризуем несколько бутстрэп-подвыборок и смотрим, насколько согласуются разбиения (по ARI). Устойчивое k — то, при котором результат воспроизводится. Это самый честный, но и самый дорогой критерий.

import numpy as np
from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler

def choose_k(X, k_range=range(2, 11), seed=0):
    """Сводная таблица диагностик по k: инерция + средний силуэт."""
    X_scaled = StandardScaler().fit_transform(X)
    rows = []
    for k in k_range:
        km = KMeans(n_clusters=k, n_init=10, random_state=seed).fit(X_scaled)
        sil = silhouette_score(X_scaled, km.labels_, sample_size=min(10_000, len(X)),
                               random_state=seed)
        rows.append({"k": k, "inertia": km.inertia_, "silhouette": round(sil, 3)})
    return rows

# Правило: сначала смотрим диагностики, потом — интерпретируемость.
# Кластеризация на 12 сегментов, которую маркетинг не может описать словами,
# бесполезна, даже если у неё лучший силуэт.

Практическая рекомендация: выбирайте k не по максимуму метрики, а по совпадению статистики и смысла. Если силуэт для k = 7 равен 0.41, а для k = 4 — 0.39, но четыре сегмента получают понятные названия («экономные новички», «лояльное ядро», …), берите четыре.

Иерархическая кластеризация: дерево вместо плоского разбиения

Иерархический подход не требует задавать k заранее. Он строит дендрограмму — дерево вложенных объединений, где по оси высоты отложено расстояние слияния. Плоское разбиение получается горизонтальным разрезом на выбранной высоте; один прогон даёт всё семейство разбиений от n кластеров до одного.

Агломеративный (снизу вверх) вариант:

AGGLOMERATIVE(X, linkage):
    каждая точка — свой кластер
    вычислить матрицу попарных расстояний D            # O(n²·d)
    повторять n−1 раз:
        (a, b) ← argmin D[a][b]
        слить a и b в новый кластер c
        пересчитать D[c][*] по формуле Ланса–Уильямса   # O(n)
    вернуть дендрограмму

Критерии связи (linkage) — здесь всё и решается

Расстояние между кластерами определяется по-разному, и от выбора зависит форма результата сильнее, чем от чего-либо ещё:

Все они — частные случаи формулы Ланса–Уильямса, позволяющей пересчитывать расстояния инкрементально, не обращаясь к исходным точкам:

d(a∪b, c) = αₐ·d(a,c) + α_b·d(b,c) + β·d(a,b) + γ·|d(a,c) − d(b,c)|

Подставляя разные (αₐ, α_b, β, γ), получаем single (½, ½, 0, −½), complete (½, ½, 0, +½), average и Ward.

Сложность

• Наивная реализация: O(n³) время, O(n²) память.
• С кучей приоритетов: O(n² log n).
• SLINK для single-linkage и CLINK для complete: O(n²) время, O(n) память (Sibson, 1973).

Квадратичная память — реальный потолок: n = 100 000 даёт матрицу расстояний на ~40 ГБ в float32. Поэтому иерархическая кластеризация в чистом виде применима к десяткам тысяч объектов; для больших данных её запускают поверх результата k-means (сначала 1000 микрокластеров, потом иерархия над их центрами) — приём, который в промышленных пайплайнах встречается постоянно.

import numpy as np
from scipy.cluster.hierarchy import dendrogram, fcluster, linkage
from scipy.spatial.distance import pdist
import matplotlib.pyplot as plt

X = np.random.default_rng(0).normal(size=(60, 4))

# Ward требует евклидовой метрики и «сырой» матрицы наблюдений
Z = linkage(X, method="ward")          # Z: (n-1, 4) — [левый, правый, высота, размер]

# Плоское разбиение: либо по числу кластеров, либо по порогу высоты
labels_k = fcluster(Z, t=4, criterion="maxclust")
labels_h = fcluster(Z, t=6.0, criterion="distance")

# Кофенетическая корреляция — насколько дендрограмма верна исходным расстояниям
from scipy.cluster.hierarchy import cophenet
c, _ = cophenet(Z, pdist(X))
print(f"кофенетическая корреляция: {c:.3f}")   # > 0.75 — дерево неплохо описывает данные

plt.figure(figsize=(11, 4))
dendrogram(Z, truncate_mode="lastp", p=20, show_leaf_counts=True)
plt.ylabel("расстояние слияния (Ward)")
plt.tight_layout()

Как читать дендрограмму. Длинная вертикальная «ножка» перед слиянием означает, что объединяемые группы были далеки друг от друга — это хорошее место для разреза. Если все слияния происходят на близких высотах, выраженной кластерной структуры в данных нет. Кофенетическая корреляция формализует это ощущение: она сравнивает расстояния в дендрограмме с исходными.

DBSCAN: кластер как область высокой плотности

DBSCAN (Density-Based Spatial Clustering of Applications with Noise, Ester, Kriegel, Sander, Xu, KDD 1996) меняет саму постановку. Кластер — это не «облако вокруг центра», а связная область, где точек на единицу объёма много, отделённая от других таких областей разрежением. Отсюда три следствия, которых нет у k-means: произвольная форма кластеров, автоматическое определение их числа и явная категория «шум».

Определения

Задаются два параметра: радиус ε и minPts.

• N_ε(p) — множество точек в шаре радиуса ε вокруг p.
• Корневая (core) точка: |N_ε(p)| ≥ minPts (обычно включая саму p).
• Граничная (border): не корневая, но попадает в ε-окрестность какой-то корневой.
• Шум (noise): ни то, ни другое.
• q прямо плотностно-достижима из p, если p — корневая и q ∈ N_ε(p). Плотностно-достижима — если есть цепочка таких шагов. Кластер = максимальное множество взаимно плотностно-связанных точек.

Именно транзитивность по корневым точкам позволяет кластеру изгибаться как угодно: он «растекается» по области, пока плотность держится выше порога.

Жизненный цикл точки при обходе

Обратите внимание на переход «шум → граничная»: точка, объявленная шумом на раннем этапе, может быть подобрана позже. А вот граничная точка, попадающая в окрестность двух разных корневых из разных кластеров, достанется тому, кто обработал её первым — DBSCAN не полностью детерминирован по граничным точкам (корневые точки и множество шума при этом определены однозначно).

Псевдокод и сложность

DBSCAN(X, eps, minPts):
    label[*] ← UNDEFINED; C ← 0
    для каждой точки p из X:
        если label[p] ≠ UNDEFINED: continue
        N ← RANGEQUERY(X, p, eps)
        если |N| < minPts:
            label[p] ← NOISE; continue
        C ← C + 1; label[p] ← C
        S ← N \ {p}                       # очередь-«семя»
        пока S не пуста:
            q ← извлечь из S
            если label[q] = NOISE: label[q] ← C      # граничная
            если label[q] ≠ UNDEFINED: continue
            label[q] ← C
            M ← RANGEQUERY(X, q, eps)
            если |M| ≥ minPts: S ← S ∪ M             # q корневая — расширяемся
    вернуть label

Всё упирается в RANGEQUERY. С пространственным индексом (KD-дерево, ball-tree, R-дерево) один запрос обходится в O(log n) при небольшом d, и общая сложность — O(n log n). Без индекса или при d ≳ 20, где деревья вырождаются, получается O(n²). Память — O(n), если не материализовать матрицу расстояний (в sklearn это управляется algorithm и n_jobs; metric="precomputed" с плотной матрицей даст O(n²) памяти).

Как подбирать eps и minPts

import numpy as np
import matplotlib.pyplot as plt
from sklearn.cluster import DBSCAN
from sklearn.neighbors import NearestNeighbors
from sklearn.preprocessing import StandardScaler

X = StandardScaler().fit_transform(raw_X)

# Эвристика: minPts >= d + 1, на практике 2*d; для шумных данных — больше.
min_pts = 2 * X.shape[1]

# k-dist график: сортированные расстояния до (minPts-1)-го соседа.
# Излом («колено») кривой — разумное eps: слева плотные области, справа шум.
nn = NearestNeighbors(n_neighbors=min_pts).fit(X)
dists, _ = nn.kneighbors(X)
kdist = np.sort(dists[:, -1])

plt.plot(kdist)
plt.xlabel("точки, отсортированные по k-dist")
plt.ylabel(f"расстояние до {min_pts}-го соседа")
plt.axhline(0.35, ls="--")     # кандидат в eps — по излому

db = DBSCAN(eps=0.35, min_samples=min_pts, n_jobs=-1).fit(X)
labels = db.labels_
n_clusters = len(set(labels)) - (1 if -1 in labels else 0)
noise_ratio = (labels == -1).mean()
print(f"кластеров: {n_clusters}, доля шума: {noise_ratio:.1%}")

Ориентиры при настройке:

• Доля шума 30–60% — ε слишком мал (или данные действительно разрежены).
• Один гигантский кластер, поглотивший всё — ε слишком велик.
• minPts управляет консервативностью: больше minPts → больше шума, чётче кластеры. Авторы рекомендуют minPts ≥ d + 1; распространённая практика — 2·d.

Главная слабость: переменная плотность

DBSCAN использует один глобальный порог плотности. Если в данных есть плотный кластер и рядом разреженный, никакое ε не устроит обоих: малое ε разорвёт разреженный на шум, большое — сольёт плотные.

Решения:

• OPTICS (Ankerst et al., 1999) — не строит разбиение, а упорядочивает точки и выдаёт reachability-график, из «долин» которого извлекаются кластеры разной плотности.
• HDBSCAN (Campello, Moulavi, Sander, 2013) — строит иерархию по всем ε сразу и выбирает наиболее устойчивые кластеры. Из параметров остаётся практически один — min_cluster_size, интуитивно понятный. С версии 1.3 доступен как sklearn.cluster.HDBSCAN. Для большинства новых задач это дефолтный выбор вместо DBSCAN.

GMM: вероятностный взгляд

k-means говорит «эта точка принадлежит кластеру 3». Смесь гауссиан (Gaussian Mixture Model) говорит «с вероятностью 0.72 — кластеру 3, с вероятностью 0.26 — кластеру 1». Это не косметика: мягкие принадлежности напрямую используются в проде для порогов уверенности и обнаружения аномалий.

Модель

Предполагаем, что данные порождены смесью K нормальных распределений:

p(x) = Σₖ πₖ · N(x | μₖ, Σₖ),   πₖ ≥ 0,  Σₖ πₖ = 1

Параметры θ = {πₖ, μₖ, Σₖ} подбираются максимизацией логарифма правдоподобия Σᵢ log p(xᵢ). Прямая максимизация невозможна аналитически из-за логарифма от суммы, поэтому вводится скрытая переменная zᵢ ∈ {1..K} — «какая компонента породила xᵢ» — и применяется EM-алгоритм (Dempster, Laird, Rubin, 1977).

EM-алгоритм

E-шаг (expectation) — считаем «ответственность» компоненты k за точку i:

γᵢₖ = πₖ N(xᵢ | μₖ, Σₖ) / Σⱼ πⱼ N(xᵢ | μⱼ, Σⱼ)

M-шаг (maximization) — обновляем параметры взвешенными оценками, где вес точки в компоненте — её ответственность:

Nₖ = Σᵢ γᵢₖ
πₖ = Nₖ / n
μₖ = (1/Nₖ) Σᵢ γᵢₖ xᵢ
Σₖ = (1/Nₖ) Σᵢ γᵢₖ (xᵢ − μₖ)(xᵢ − μₖ)ᵀ

Гарантируется, что log-правдоподобие не убывает на каждой итерации (это следует из неравенства Йенсена и свойств нижней ELBO-границы), но, как и у k-means, сходимость — к локальному максимуму.

Связь с k-means. Зафиксируйте Σₖ = σ²I и устремите σ → 0. Ответственности вырождаются в 0/1 — вся масса уходит ближайшему центру, и EM превращается ровно в алгоритм Ллойда. k-means — это предельный частный случай GMM с жёстким назначением и сферическими равными ковариациями. Отсюда же понятно, что даёт GMM сверх k-means: полная ковариационная матрица позволяет кластерам быть вытянутыми и повёрнутыми.

Практика

import numpy as np
from sklearn.mixture import GaussianMixture

# covariance_type определяет гибкость и число параметров на компоненту:
#   'spherical' — K*(d+2)      : шары разного радиуса (≈ k-means)
#   'diag'      — K*(2d+1)     : оси эллипсоидов параллельны координатным
#   'tied'      — общая Σ      : одинаковая форма у всех компонент
#   'full'      — K*(d(d+3)/2+1): произвольные эллипсоиды (по умолчанию)
gm = GaussianMixture(
    n_components=4,
    covariance_type="full",
    n_init=10,               # EM тоже застревает в локальных максимумах
    reg_covar=1e-6,          # защита от вырождения ковариации
    random_state=0,
).fit(X)

hard = gm.predict(X)              # жёсткие метки
soft = gm.predict_proba(X)        # (n, K) — ответственности
logp = gm.score_samples(X)        # log p(x): низкие значения = аномалии

# Выбор числа компонент по информационным критериям.
# BIC штрафует сложность сильнее AIC и обычно даёт более экономную модель.
scores = []
for k in range(1, 11):
    m = GaussianMixture(n_components=k, covariance_type="full",
                        n_init=5, random_state=0).fit(X)
    scores.append((k, m.bic(X), m.aic(X)))
best_k = min(scores, key=lambda t: t[1])[0]

Сложность одной итерации EM: O(n·K·d²) для full (из-за обращения ковариаций — O(K·d³) разово через разложение Холецкого) и O(n·K·d) для diag. На высокой размерности full быстро становится неподъёмным и переобучается: d = 100 и K = 10 — это уже 50 тысяч параметров ковариаций. Разумный порядок действий: PCA до 10–50 компонент, потом GMM.

Вырождение — специфическая болезнь GMM: компонента «садится» на одну точку, её ковариация стремится к нулю, а правдоподобие — к бесконечности. Формально это глобальный максимум, практически — мусор. Защита: reg_covar (добавка εI к диагонали) и байесовский вариант BayesianGaussianMixture с априорами Дирихле, который вдобавок умеет сам «выключать» лишние компоненты, обнуляя их веса.

Оценка качества без разметки

Первое, что нужно осознать: accuracy к кластеризации неприменима. Метки кластеров — произвольные идентификаторы; перенумерация не меняет разбиение, но обнуляет accuracy. Нужны метрики, инвариантные к перестановке меток.

Внутренние (только по X и меткам)

Общая оговорка: все внутренние метрики построены на предположениях о форме кластеров и потому систематически занижают оценку DBSCAN/HDBSCAN на данных сложной формы. Сравнивать по силуэту k-means и DBSCAN на полумесяцах — методологическая ошибка: силуэт выберет геометрически неправильное решение.

Внешние (когда есть эталонная разметка)

from sklearn.metrics import (adjusted_rand_score, adjusted_mutual_info_score,
                             silhouette_score, calinski_harabasz_score,
                             davies_bouldin_score)

def report(X, labels, y_true=None):
    """Сводка качества. mask нужен, чтобы шум DBSCAN (-1) не портил внутренние метрики."""
    out = {}
    mask = labels != -1
    if len(set(labels[mask])) > 1:
        out["silhouette"] = silhouette_score(X[mask], labels[mask])
        out["calinski_harabasz"] = calinski_harabasz_score(X[mask], labels[mask])
        out["davies_bouldin"] = davies_bouldin_score(X[mask], labels[mask])
    out["noise_ratio"] = float((~mask).mean())
    if y_true is not None:
        out["ARI"] = adjusted_rand_score(y_true, labels)   # инвариантна к перестановке меток
        out["AMI"] = adjusted_mutual_info_score(y_true, labels)
    return out

Тонкость с шумом: считать силуэт, включая точки -1, некорректно — «шум» это не кластер. Но и оценивать только по кластеризованной части нечестно: алгоритм, объявивший шумом 70% данных, легко получит отличный силуэт на оставшихся 30%. Всегда сообщайте noise_ratio рядом с метрикой.

Самая надёжная оценка — внешняя проверка полезности. Сегменты пользователей проверяются A/B-тестом на разных коммуникациях; кластеры логов — тем, находит ли по ним дежурный инженер реальные инциденты. Внутренние метрики — только для отсева заведомо плохих конфигураций.

Сравнение: что выбирать

Сводная таблица для принятия решения:

Последняя строка часто оказывается решающей в проде. k-means и GMM дают обученную модель, которую можно сериализовать и применять к потоку новых объектов за O(k·d). DBSCAN и иерархическая кластеризация — это разбиение конкретной выборки, а не модель; чтобы разметить новый объект, нужен либо пересчёт, либо надстройка (обучить классификатор на полученных метках — распространённый и вполне легитимный приём).

Типичные ошибки

1. Не масштабировали признаки перед k-means/DBSCAN/иерархической. Результат определяется признаком с наибольшим разбросом.
2. Кластеризация «в лоб» на сотнях признаков. В высокой размерности расстояния концентрируются; сначала PCA/UMAP, потом кластеризация.
3. One-hot от категории с 500 уровнями + евклидово расстояние. Расстояние между любыми двумя разными категориями одинаково (√2) — информации почти нет. Используйте k-modes, Gower или целевое/частотное кодирование.
4. Интерпретация k-means-кластеров как «настоящих групп». k-means всегда вернёт k кластеров, даже на равномерном шуме. Проверяйте gap statistic или тест Хопкинса на наличие структуры вообще.
5. Выбор k только по локтю. На реальных данных излом обычно размыт; смотрите несколько диагностик плюс интерпретируемость.
6. Сравнение алгоритмов по силуэту. Силуэт встроенно предпочитает сферические кластеры и потому нечестен к плотностным методам.
7. Игнорирование доли шума в DBSCAN. 60% шума — это не «нашли чистые кластеры», а неправильный ε.
8. Одна кластеризация навсегда. Данные дрейфуют; сегменты полугодовой давности могут описывать несуществующую реальность. Нужен мониторинг стабильности — см. MLOps.
9. Забыли, что DBSCAN недетерминирован по граничным точкам — и удивляются, почему метки «прыгают» между запусками при разном порядке строк.
10. Кластеризация вместо классификации. Если разметка есть или её можно получить, обучение с учителем почти всегда точнее. Кластеризация — для случая, когда классов никто не знает.

Как это применяют в продакшене

Сегментация клиентов. Классика — RFM (Recency, Frequency, Monetary): три признака, логарифмирование (распределения тяжелохвостые), стандартизация, k-means на k = 4…8. Результат отдаётся маркетингу как набор именованных сегментов, каждый со своей коммуникацией. Критерий успеха — не силуэт, а разница конверсии в A/B-тесте.

Аномалии и мониторинг. DBSCAN/HDBSCAN на признаках сессии: всё, что помечено -1, — кандидаты в фрод или сбои. GMM даёт непрерывную альтернативу: score_samples возвращает log p(x), и нижний перцентиль — аномалии. Плюс подхода в том, что порог настраивается под допустимую нагрузку на дежурную смену.

Векторное квантование и ANN-поиск. Индексы IVF в FAISS строятся именно k-means: пространство разбивается на ячейки Вороного, поиск идёт только по нескольким ближайшим ячейкам вместо полного перебора. Product Quantization — тот же k-means, применённый к подвекторам, — сжимает 768-мерный float32-эмбеддинг (3 КБ) до 96 байт. Без этого векторные БД на миллиарды объектов не существовали бы.

Тематическое моделирование эмбеддингов. Современный стандарт (например, BERTopic) — конвейер «эмбеддинги предложений → UMAP до 5–10 измерений → HDBSCAN → c-TF-IDF для описания тем». HDBSCAN здесь важен тем, что не обязан приписывать тему каждому документу: мусор честно уходит в шум.

Дедупликация и entity resolution. Блокирование кандидатов (LSH/MinHash) → граф похожести → связные компоненты или иерархическая кластеризация с порогом. Single-linkage тут одновременно нужен (транзитивность дублей) и опасен (chaining склеит два разных объекта через одну ошибочную пару) — поэтому порог подбирают на размеченном наборе пар. Родственный приём — полуразметка: кластеризуем непомеченный корпус, размечаем по несколько объектов из каждого кластера и распространяем метки, получая стартовый датасет для supervised-модели.

Инженерная сторона. Кластеризация в проде — почти всегда батч-джоб (ежедневный/еженедельный) с сохранением модели и версии сегментов. Типичные грабли: перенумерация кластеров между запусками (кластер «1» вчера и сегодня — разные группы). Лечение — сопоставление новых кластеров со старыми по расстоянию между центроидами и стабильные бизнес-идентификаторы сегментов.

Мини-итог

• Кластеризация не находит «истинную» структуру — она находит структуру, соответствующую своему определению кластера. Теорема Клейнберга объясняет, почему универсального алгоритма нет.
• k-means минимизирует WCSS алгоритмом Ллойда, требует k, даёт выпуклые кластеры, стоит O(nkdi) — базовый выбор для больших числовых данных при разумной геометрии.
• Иерархическая строит дендрограмму и не требует k, но O(n²) по памяти; выбор linkage важнее выбора всего остального.
• DBSCAN определяет кластер через плотность: произвольные формы, автоматическое k, явный шум; страдает от переменной плотности — используйте HDBSCAN.
• GMM — вероятностная модель с мягкими принадлежностями и эллипсоидальными кластерами; k-means является её вырожденным случаем; число компонент выбирается по BIC.
• Масштабирование обязательно, высокая размерность губительна, внутренние метрики смещены в пользу сферических кластеров, а окончательный судья — полезность сегментов в задаче.

Источники

• Ester, Kriegel, Sander, Xu. A Density-Based Algorithm for Discovering Clusters, KDD 1996 — оригинальная статья DBSCAN.
• Arthur, Vassilvitskii. k-means++: The Advantages of Careful Seeding, SODA 2007.
• Kleinberg. An Impossibility Theorem for Clustering, NIPS 2002.
• Tibshirani, Walther, Hastie. Estimating the number of clusters via the gap statistic, JRSS-B 2001.
• Rousseeuw. Silhouettes: a graphical aid to the interpretation and validation of cluster analysis, 1987.
• Campello, Moulavi, Sander. Density-Based Clustering Based on Hierarchical Density Estimates, PAKDD 2013 — HDBSCAN.
• Hastie, Tibshirani, Friedman. The Elements of Statistical Learning, глава 14 — бесплатный PDF.
• Bishop. Pattern Recognition and Machine Learning, глава 9 — самый аккуратный разбор EM и GMM.
• scikit-learn: Clustering — включая знаменитую сравнительную картинку поведения алгоритмов на игрушечных наборах.
• SciPy: hierarchical clustering и документация HDBSCAN.

Смежные материалы трека: расстояния и метрические методы разбираются в статье про k ближайших соседей; подготовка и масштабирование признаков — в данных и признаках; ковариационные матрицы и разложения, лежащие в основе GMM и PCA, — в математике для ML; метрики и валидация — в оценке моделей.

Что дальше

Мы дважды упирались в одно и то же ограничение: в высокой размерности расстояния перестают различать объекты, а значит все рассмотренные алгоритмы деградируют. Кроме того, кластеры удобно смотреть глазами, а десятимерное пространство глазами не смотрится. Оба вопроса решает следующая статья: Снижение размерности: PCA, SVD, t-SNE, UMAP.