DBSCAN

За 2015–2025 годы на Земле случилось 19 073 землетрясения магнитудой от 5. На карте видно, что они не рассыпаны как попало: они тянутся цепочками вдоль границ литосферных плит, сгущаются в рои и изредка случаются посреди плиты в одиночку. Сколько тут групп? k-средних нужно назвать число заранее, и он нарежет точки на круглые куски. DBSCAN число групп не спрашивает. Он ищет плотные места, где у точки много соседей, соединяет их цепочками в кластеры любой формы, а одинокие толчки оставляет в стороне, серым шумом. При радиусе 100 км и пороге в десять точек получается 162 кластера, и 14 % землетрясений остаются шумом.

Подвигайте радиус. При 25 км больше половины землетрясений становятся шумом, и остаются только самые плотные рои. При 400 км кластеры сливаются в несколько гигантских поясов. Переключатель «и не дальше 30 дней» делает соседями только те толчки, которые были близко и по месту, и по времени. Тогда вместо границ плит видны рои: серии повторных толчков вокруг сильных землетрясений. Кнопка k-средних раскрашивает те же точки на столько же групп по ближайшему центру, и разница с DBSCAN видна сразу: пояса порезаны на куски примерно одного размаха, а у одиноких толчков посреди океана тоже появляется группа.

Идея

Кластер — это место, где точек густо, а между кластерами — пусто. Густоту меряют самым простым способом: сколько точек лежит в круге радиуса ε вокруг точки. Если много — точка внутри кластера, и её соседи тоже в кластере. Если мало — она на краю кластера или вовсе одна.

Параметров два: радиус ε и порог minPts. Соседи точки p — все точки на расстоянии не больше ε от неё, включая её саму. Обозначим их N(p). Каждая точка получает одну из трёх ролей.

• Ядро — точка, у которой |N(p)| ≥ minPts. Ядра — плотная внутренность кластеров.
• Край — не ядро, но в радиусе ε от какого-нибудь ядра. Край входит в кластер, но кластер через него дальше не растёт.
• Шум — всё остальное: точка не ядро, и ни одного ядра рядом нет.

Кластер — это ядра, связанные цепочками, в которых каждые два соседних ядра — соседи, вместе с краями этих ядер. Собирает их алгоритм так. Берём ещё не просмотренную точку и считаем её соседей. Если их меньше minPts, точка пока считается шумом: позже её может подобрать ядро, и тогда она станет краем. Если соседей хватает, начинается новый кластер: все соседи встают в очередь. Каждую точку из очереди отправляем в этот кластер, и если она тоже ядро, в очередь встают и её соседи. Когда очередь опустеет, кластер готов, и мы переходим к следующей непросмотренной точке. По сути это обход в ширину по графу, где рёбра идут от ядер к их соседям.

На теореме держится и задача о кластерах на карте эпицентров: число кластеров — это число компонент связности в графе «ядро — соседнее ядро», и его можно собрать любым способом, хоть системой непересекающихся множеств. Если же minPts = 2, ядром становится любая точка, у которой есть хоть один сосед. Тогда кластеры DBSCAN — те же куски, на которые распадается сеть рёбер короче ε, то есть одиночная связь из статьи о минимальном остовном дереве, а точки без соседей уходят в шум. Порог minPts больше двух защищает от цепочек: тонкий мостик из редких точек не соединяет два кластера, если точки мостика не ядра.

По шагам

На целом мире шагов не разглядеть, поэтому возьмём одну дугу островов на северо-западе Тихого океана: 185 землетрясений за 2018–2021 годы на плоской карте этой области. DBSCAN идёт по ним в порядке времени. Круг вокруг текущей точки — её радиус ε. Ядра — большие кружки, края — кружки с пустой серединой, шум — серые точки.

При ε = 100 км и minPts = 4 получается 6 кластеров: два длинных, по 67 и 59 точек, вытянуты вдоль дуги, ещё один из 35 точек, и три маленьких, по 4 точки. Из 185 точек 167 — ядра, 6 — края и 12 — шум. Увеличьте ε до 120 км, и северный кластер из 67 точек сольётся со средним из 35 и одной из маленьких групп в кластер из 107 точек: при большем радиусе между ними нашлась цепочка ядер. Поднимите minPts до 5, и маленькие группы рассыплются в шум, а кластеров останется три. Переключатель «k-средних» раскрашивает те же точки на столько же групп по ближайшему центру. Длинную дугу он режет на куски примерно одного размаха, шума у него нет вовсе, а одну из шести групп он собирает из разрозненных одиночек, которых DBSCAN счёл шумом.

Код

Землетрясения лежат на сфере, и расстояние между ними — дуга по поверхности Земли. Широта и долгота для расстояний неудобны: на шестидесятой широте градус долготы вдвое короче, чем на экваторе, а на 180-м меридиане долгота прыгает со 180 на −180. Поэтому каждая точка превращается в единичный вектор (x, y, z), как в статье о k-d дереве. Дуге d на Земле радиуса R соответствует хорда 2·sin(d/2R), и хорда растёт вместе с дугой. Значит, «дуга не больше ε» — то же, что «хорда не больше хорды ε», и сравнивать можно обычные расстояния в трёхмерном пространстве.

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

19073 earthquakes of magnitude 5+, 2015-2025
100 km, 10 points: 162 clusters, 15224 cores, 1177 border points, 2672 noise (14.0%); largest 2917, 2473, 2325; 2373291 distances instead of 181880128
100 km and 30 days, 10 points: 126 clusters, 3499 cores, 362 border points, 15212 noise (79.8%); largest 465, 313, 215; 2373291 distances instead of 181880128

Вместо 181,9 млн пар сетка посчитала 2,4 млн расстояний — в 77 раз меньше. Три самых больших кластера — пояса вдоль глубоководных желобов, где одна плита уходит под другую. Самый большой, 2917 землетрясений, тянется на 4779 км по островам на юго-западе Тихого океана. Второй, 2473 землетрясения, — дуга на 7469 км по северо-западу и северу Тихого океана: несколько островных дуг подряд, сцепленных ядрами в одну цепь. Третий, 2325, — узкий жёлоб длиной около 2500 км к югу от экватора, по обе стороны от 180-го меридиана.

Соседи во времени

Вторая строка вывода — тот же алгоритм с другим понятием соседа: две точки соседи, если между ними не больше 100 км и не больше 30 дней. DBSCAN не важно, что такое расстояние. Ему нужен только ответ на вопрос «соседи ли эти две точки», и любое правило, которое на него отвечает, даёт свою кластеризацию. Здесь кластеры — уже не границы плит, а рои: сильное землетрясение и серия повторных толчков после него, плотная и в пространстве, и во времени. Шумом становятся 79,8 % точек: у большинства сильных землетрясений нет и девяти соседей ближе 100 км и 30 дней. Самый большой рой — 465 толчков магнитудой от 5 за пять месяцев 2025 года на северо-западе Тихого океана; самый сильный из них — магнитудой 8,8.

Сетка при этом не меняется: соседи по времени — подмножество соседей по месту, поэтому кубики по-прежнему отсекают далёкие точки, а время проверяется вторым условием. Число сравнений то же, 2 373 291.

Как выбрать ε и minPts

Порог minPts выбрать проще: он сглаживает оценку плотности, и в шумных данных его поднимают. Авторы DBSCAN предлагали для точек на плоскости брать 4. В 1998 году Йорг Зандер с соавторами предложили удвоенную размерность данных, 2·d, а для больших и шумных наборов советуют брать ещё больше. С ε труднее. Если есть знание о предмете, его и берут. Эрих Шуберт с соавторами в статье 2017 года приводят такой пример: для координат GPS специалист может решить, что соседи — точки ближе километра. Наши 100 км — просто круглое число, и график ниже показывает, что оно значит для этих данных.

Для каждой точки найдём расстояние до её (minPts − 1)-го ближайшего соседа — это её расстояние ядра. Минус один — потому что сама точка уже первая из minPts. Точка — ядро ровно при тех ε, которые не меньше этого расстояния. Выстроим точки по убыванию расстояния ядра: получится отсортированный график расстояний до соседей, предложенный ещё в исходной статье о DBSCAN. Горизонтальная линия на высоте ε делит его на две части. Всё, что ниже линии, — ядра, всё, что выше, — края или шум.

Авторы DBSCAN советовали искать на этом графике первую «впадину», где кривая резко переламывается, а если известно, какая доля данных — шум, отложить её слева и взять ε в этом месте. На настоящих данных излома часто нет. Здесь его тоже нет: кривая на логарифмической шкале спускается плавно. При minPts = 10 у половины землетрясений девятый сосед не дальше 40 км, а у тысячи самых одиноких — дальше 225 км. При ε = 100 км ядрами становятся 15 224 точки, 79,8 %, и график в этом месте говорит то же, что листинг. Ещё 1177 точек оказываются краями, так что в шум уходят 14 %.

В той же статье Шуберт с соавторами называют два признака неудачных параметров. Первый — доля шума: разумной обычно бывает доля от 1 до 30 %. Второй — размер самого большого кластера: если в нём больше 20–50 % всех точек, попавших в кластеры, ε, скорее всего, слишком велик. На землетрясениях при 100 км шум — 14 %, а самый большой кластер — 17,8 % точек в кластерах. При 150 км самый большой кластер — уже 35 %, при 300 км — 64 %: пояса сливаются. При 50 км шум — 31,6 %, чуть выше верхней границы. Оба признака сужают выбор примерно до 50–150 км. Там же оговорена мелочь, на которой часто ошибаются. Соседи в поиске ближайших соседей не включают саму точку, а в подсчёте плотности DBSCAN — включают. Поэтому для порога minPts нужен (minPts − 1)-й сосед. Этот сдвиг на единицу — подвох задачи о пороге плотности.

DBSCAN и k-средних

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

В сцене с картой k-средних получает то же число групп, 162, и стартует с k-means++. Поперечник типичной группы k-средних (медиана) — 1114 км, типичного кластера DBSCAN — 211 км. Зато самый длинный кластер DBSCAN тянется на 7469 км, а самая длинная группа k-средних — на 2542 км: пояс, который DBSCAN видит целиком, k-средних режет на куски. 19 кластеров DBSCAN он делит на несколько заметных частей, а 52 его группы смешивают точки разных кластеров DBSCAN. И у k-средних нет шума: все 2672 одиноких толчка он раздал ближайшим группам. До центра своей группы у них по медиане 375 км, а у 64 — больше тысячи.

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

Сколько это стоит

DBSCAN спрашивает соседей каждой точки ровно один раз: точка уходит из «непросмотренных», как только её соседей посчитали. Поэтому время работы — n запросов соседей плюс время, чтобы перебрать всех найденных соседей: O(n·Q + Σ|N(p)|), где Q — цена одного запроса. Без индекса запрос — перебор всех точек, Q = O(n), и всё стоит O(n²). С сеткой запрос смотрит только 27 кубиков, и на землетрясениях это в 77 раз меньше расстояний. Вместо сетки годится и k-d дерево с поиском всех точек в радиусе. В scikit-learn у DBSCAN есть выбор: k-d дерево, шаровое дерево или перебор.

Перебор пар растёт как n², как и положено. Сетка в Python быстрее перебора в десятки раз, но и её линия идёт круче прямой пропорции. Здесь нет ошибки. Землетрясения взяты по порядку времени, и за вдвое больший срок в тех же 100 км накапливается вдвое больше соседей. Слагаемое Σ|N(p)| растёт как квадрат, и никакой индекс его не уберёт: найденных соседей всё равно нужно перебрать. Сетка убирает только сравнения с далёкими точками. В Go разрыв меньше, в 5–16 раз: там перебор пар — тугой цикл, а сетка платит за поиск кубиков в словаре.

Худший случай остаётся квадратичным при любом индексе: если все точки лежат в радиусе ε друг от друга, каждый запрос возвращает почти всех. В исходной статье 1996 года время DBSCAN с R*-деревом оценено как O(n log n) в среднем, в расчёте на то, что окрестности маленькие. В 2013 году Аде Гунаван в магистерской работе заметил, что в худшем случае алгоритм квадратичный, и для плоскости построил точный алгоритм за O(n log n) на сетке с клетками ε/√2. Любые две точки такой клетки — соседи, и если в клетке набралось minPts точек, все они ядра без единого сравнения. В 2015 году Цзюньхао Гань и Юфэй Тао показали, что в трёх и более измерениях точный DBSCAN, по-видимому, требует порядка n^4/3 времени. Обойти эту оценку можно, только решив быстрее одну из давних трудных задач вычислительной геометрии, а в это мало кто верит. Взамен они предложили приближённый DBSCAN за O(n) в среднем. Через два года Шуберт, Зандер, Эстер, Кригель и Сюй ответили: упрекать стоило не DBSCAN, а надежду на пространственные индексы. На одном из наборов данных из той работы при наименьшем ε, взятом в опытах Ганя и Тао, больше 99 % точек уже лежат в одном кластере. При разумно выбранных ε обычный DBSCAN с хорошим индексом работает не хуже нового приближённого метода.

Где алгоритм ошибается

• Одна плотность на всех. Радиус ε один на все данные. Если рядом есть густой рой и разреженный пояс, ε, подобранный для пояса, сольёт рой с соседями, а подобранный для роя превратит пояс в шум. Это видно на карте: при 50 км самые плотные жёлоба ещё держатся целыми кластерами, а срединные хребты посреди океанов почти целиком ушли в шум. Эту слабость лечат OPTICS (1999), который одним проходом упорядочивает точки так, что из порядка читаются кластеры для всех ε до заданного предела, и HDBSCAN (2013), который выбирает кластеры из целой иерархии плотностей.
• Цепочки ядер. Если два плотных места соединяет вереница ядер, они станут одним кластером, хотя на глаз это два. Так собран второй по величине кластер землетрясений, дуга на 7469 км. minPts защищает от мостиков из редких точек, но не от мостика из ядер.
• Края зависят от порядка. Край, у которого рядом ядра двух кластеров, достаётся тому, кто дотянулся первым. Поэтому кластер может оказаться даже меньше minPts: часть соседей его первого ядра уже забрал кластер, собранный раньше. В PostGIS документация функции ST_ClusterDBSCAN прямо предупреждает, что такие края распределяются произвольно, и советует задать порядок обхода.
• Много измерений. В пространстве с десятками признаков расстояния между любыми точками становятся похожими, и «плотное место» теряет смысл. Перед DBSCAN размерность часто снижают, например методом главных компонент, а сетка в высоких размерностях бесполезна: соседних клеток 3^d.
• Масштаб признаков. Расстояние смешивает все координаты, и признак в больших единицах заглушит остальные. В соседях «по месту и времени» эта проблема решена двумя отдельными порогами — 100 км и 30 дней, — а не одной формулой вроде «километры плюс дни».

История

DBSCAN предложили Мартин Эстер, Ханс-Петер Кригель, Йорг Зандер и Сяовэй Сюй на конференции KDD в 1996 году. Статья называлась «Алгоритм кластеризации по плотности для поиска кластеров в больших пространственных базах данных с шумом», и проверяли алгоритм на точках плоскости — искусственных и настоящих географических из набора данных SEQUOIA 2000. В 1998 году Зандер, Эстер, Кригель и Сюй обобщили его до GDBSCAN: соседство в нём задаётся любым симметричным правилом «соседи или нет», а кластеризовать можно не только точки, но и протяжённые объекты, например многоугольники на карте. В 1999 году Анкерст, Бройниг, Кригель и Зандер описали OPTICS, в 2013 году Кампелло, Мулави и Зандер — HDBSCAN. В 2014 году исходная статья получила награду SIGKDD «Испытание временем» за влияние спустя годы. DBSCAN есть в scikit-learn, ELKI, R и Weka, а с версии 2.3 — и в PostGIS, прямо в запросах к базе данных.

Задачи

Первая задача — сам DBSCAN на плоскости с сеткой: сколько кластеров, ядер, краёв и шума. Вторая — про выбор ε: как сразу ответить, сколько ядер будет при любом радиусе.