Кластеризация – одно из самых распространенных применений методов машинного обучения без учителя. Задача кластеризации возникает, когда вам нужно разделить данные (или объекты, описываемые каждый неким набором параметров/координат) на группы объектов, похожих/близких между собой.

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

Другой популярный метод деления кластеризаций – на плоские (которые получают одно деление на кластера) и иерархические методы (которые получают целое семейство кластеризаций). К иерархическим методам относятся в основном аггломеративные и дивизивные методы.

Каждый метод кластеризации основан на какой-то своей модели взаимного расположения объектов (хотя эта модель не всегда конкретизируется). Поэтому, применяя разные методы кластеризации к одним и тем-же данным вы будете получать различное разделение одних и тех-же объектов на кластера. Интересную (и часто приводимую картинку) можно увидеть здесь: https://scikit-learn.org/stable/modules/clustering.html.

Наиболее простыми и понятными алгоритмами являются плоские алгоритмы: смесь гауссов (Gaussian Mixture) и DBSCAN. Рассмотрим их чуть более подробно, чтобы понять почему они так часто используются, в чем их плюсы и минусы, и почему мне захотелось их объединить.

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

Для поиска значения гиперпараметра кластеризации можно применять различные методы. Можно просто выбрать его значение в соответствии со своим пониманием того, какие там объекты и сколько их типов в вашем датасете. Это называется внешним критерием или использованием эксперта. А можно выбрать значение гиперпараметра таким образом, чтобы оптимизировать (минимизировать или максимизировать) значение некоторой функции – метрики качества кластеризации. Метрик качества тоже много разных, и каждый из них опирается на какую-то модель взаимного расположения объектов. Оптимизация метрики качества кластеризации – это так-называемый внутренний критерий оценки качества.

Методы, основанные на расстоянии: DBSCAN+Silhouette

Основу метода DBSCAN составляет расстояние между объектами (точками в N-мерном пространстве). В библиотеке Sklearn существует две реализации этого алгоритма – через анализ координат объектов (функция DBSCAN) и через анализ матрицы расстояний между ними (функция dbscan).

Одной из распространенных метрик является Силуэт (Silhouette) – фактически это отношение расстояния между внешними границами кластеров к расстоянию между центрами кластеров. Чем меньше размеры кластеров – тем ближе силуэт к единице. Обычно он хорош, когда все кластера имеют форму шара (т.е. линейный размер кластера не зависит от направления, в котором он вычисляется) или когда эти размеры значительно меньше расстояний между кластерами (кластера можно заключить в непересекающиеся сферы).

Посмотрим, как работает DBSCAN+Silhouette

from sklearn.datasets import make_blobs
import matplotlib.pyplot as plt
import numpy as np
from sklearn.metrics import silhouette_score
from sklearn.cluster import DBSCAN 
S,N=1000,5
X,Y=make_blobs(n_samples=S, n_features=2, centers=N, cluster_std=0.3, random_state=42)
fig,axs=plt.subplots(1,3,figsize=(7,2.5),constrained_layout=True)
axs[0].set_title('source clusters')
axs[0].scatter(X[:,0],X[:,1],c=Y,s=1)
arr=[]
for eps in np.linspace(0.1,10,100):
 clf=DBSCAN(eps=eps)
 Yp=clf.fit_predict(X)
 ssc=silhouette_score(X,Yp)
 arr.append([eps,ssc])
arr=np.array(arr)
axs[1].set_title('Силуэт')
axs[1].set_xlabel('eps')
axs[1].set_ylabel('Silhouette')
axs[1].plot(arr[:,0],arr[:,1])
eps=arr[np.argmax(arr[:,1]),0]
clf=DBSCAN(eps=eps)
Yp=clf.fit_predict(X)
axs[2].set_title('Silh-optimal clusters')
axs[2].scatter(X[:,0],X[:,1],c=Yp,s=1);

В результате получим картинку:

Рис.1. Случай с изолированными кластерами

Слева показаны оригинальные точки датасета, цвета обозначают различные кластера, по центру - зависимость коэффициента силуэта от \eps - гиперпараметра кластеризации DBSCAN, справа - результат кластеризации с оптимальным значением \eps, соответствующим максимуму коэффициента силуэта. Видно, что кластеризация неплохо разделяет изолированные кластера.

Усложним задачу: увеличим размер кластеров (заменим в 7й строчке кода параметр cluster_std=0.3 на cluster_std=0.7) и снова запустим код:

Рис.2. Случай DBSCAN+Silhouette с соприкасающимися кластерами

Видно, что два близко расположенных (соприкасающихся) кластера перестали разделяться.

Еще эффективнее эту проблему можно продемонстрировать добавлением шума:

X,Y=make_blobs(n_samples=S, n_features=2, centers=N, cluster_std=0.3, random_state=42)
Xn=(np.random.rand(S,2)-0.5)*(X.max()+abs(X.min()))
Yn=np.array([Y.max()+1]*S)
X=np.concatenate([X,Xn])
Y=np.concatenate([Y,Yn])
fig,axs=plt.subplots(1,3,figsize=(7,2.5),constrained_layout=True)
axs[0].set_title('source clusters')
axs[0].scatter(X[:,0],X[:,1],c=Y,s=1)
arr=[]
for eps in np.linspace(0.1,10,100):
 clf=DBSCAN(eps=eps)
 Yp=clf.fit_predict(X)
 if np.unique(Yp).shape[0]>1:
  ssc=silhouette_score(X,Yp)
  arr.append([eps,ssc])
arr=np.array(arr)
axs[1].set_title('Силуэт')
axs[1].set_xlabel('eps')
axs[1].set_ylabel('Silhouette')
axs[1].plot(arr[:,0],arr[:,1])
eps=arr[np.argmax(arr[:,1]),0]
clf=DBSCAN(eps=eps)
Yp=clf.fit_predict(X)
axs[2].set_title('Silh-optimal clusters')
axs[2].scatter(X[:,0],X[:,1],c=Yp,s=1);

В этом случае кластера перестали разделяться вообще:

Методы, основанные на модели распределения: GM+BIC

Основу метода GM (Gaussian Mixture) составляет модель распределения точек в одном кластере – многомерное нормальное (т.е. гауссово) распределение. Модель предполагает, что вероятность появления точки в пространстве описывается суммой нормальных распределений, каждое из которых взято со своим весом и со своими параметрами (ширина распределения, наклон осей). Каждое такое гауссово распределение в многомерной системе координат интуитивно может быть представлено в виде эллипсоида.

Одной из распространенных метрик качества гауссовой смеси является BIC (байесовский информационный критерий) – фактически это величина обратная правдоподобию - полной вероятности верного описания моделью всех экспериментальных точек (чем вероятность выше – тем BIC меньше). Минимум этой функции будет соответствовать наилучшей модели. Очевидно, что если число параметров модели сравнимо с количеством объектов, то такой моделью можно описать объекты абсолютно точно. Поэтому в критерии BIC есть небольшая аддитивная добавка (штраф) – чем больше в модели параметров, тем больше штраф за использование такой модели. Это не позволяет использовать модели со слишком большим количеством параметров. Эту добавку можно интерпретировать, как принцип Оккама: он не дает возможности увеличивать число параметров модели без необходимости и использовать модели со слишком большим количеством параметров. Выбрать эту добавку можно различными способами, и существует даже научная статья [https://onlinelibrary.wiley.com/doi/abs/10.1111/1468-0262.00273], где показывается, что в определенных ограничениях все эти добавки приводят к одинаковому числу оптимальных кластеров.

Посмотрим, как работает GM+BIC:

from sklearn.datasets import make_blobs
import matplotlib.pyplot as plt
import numpy as np
from sklearn.mixture import GaussianMixture 
S,N=1000,5
X,Y=make_blobs(n_samples=S, n_features=2, centers=N, cluster_std=0.3, random_state=42)
fig,axs=plt.subplots(1,3,figsize=(7,2.5),constrained_layout=True)
axs[0].set_title('source clusters')
axs[0].scatter(X[:,0],X[:,1],c=Y,s=1)
arr=[]
for n_cl in range(1,20):
 clf=GaussianMixture(n_components=n_cl)
 clf.fit(X)
 bic=clf.bic(X)
 arr.append([n_cl,bic])
arr=np.array(arr)
axs[1].set_title('BIC')
axs[1].set_xlabel('n_components')
axs[1].set_ylabel('BIC')
axs[1].plot(arr[:,0],arr[:,1])
n_cl=arr[np.argmin(arr[:,1]),0]
clf=GaussianMixture(n_components=int(n_cl))
Yp=clf.fit_predict(X)
axs[2].set_title('BIC-optimal clusters')
axs[2].scatter(X[:,0],X[:,1],c=Yp,s=1)

В качестве гиперпараметра GM выступает число кластеров, и в случае простых изолированных кластеров мы получим хороший результат:

Рис.4. Случай GM+BIC с изолированными кластерами

В случае соприкасающихся кластеров или зашумленных данных он также даст неплохие результаты:

Рис.5. Случай GM+BIC с соприкасающимися и зашумленными кластерами

Если GM такой устойчивый, почему его не использовать всегда?

Потому-что он работает только с кластерами эллиптической формы. Если форма кластеров отличается от эллиптической, начинаются проблемы. Сгенерируем датасет посложнее:

def MakeCircles():
 NN=1000
 r=np.random.rand(NN)*10+50
 r2=np.random.rand(NN)*10
 r3=np.random.rand(NN)*10+80
 phi=np.random.rand(NN)*6.28
 phi2=np.random.rand(NN)*6.28
 phi3=np.random.rand(NN)*6.28
 X1=np.zeros((NN,2))
 Y1=np.zeros(NN)
 X1[:,0]=r*np.cos(phi)
 X1[:,1]=r*np.sin(phi)
 Y1[:]=1

 X2=np.zeros((NN,2))
 Y2=np.zeros(NN)
 X2[:,0]=r2*np.cos(phi2)
 X2[:,1]=r2*np.sin(phi2)
 Y2[:]=0

 X3=np.zeros((NN,2))
 Y3=np.zeros(NN)
 X3[:,0]=r3*np.cos(phi3)
 X3[:,1]=r3*np.sin(phi3)
 Y3[:]=0

 X=np.concatenate([X1,X2,X3],axis=0)
 Y=np.concatenate([Y1,Y2,Y3],axis=0)
 idx=np.array(range(X.shape[0]))
 np.random.shuffle(idx)
 X=X[idx]
 Y=Y[idx]

 Nn=300
 Xn1=np.random.rand(Nn,1)*(X[:,0].max()-X[:,0].min())+X[:,0].min()
 Xn2=np.random.rand(Nn,1)*(X[:,1].max()-X[:,1].min())+X[:,1].min()
 Xn=np.concatenate([Xn1,Xn2],axis=1)
 Yn=np.zeros(Nn)
 Yn[:]=Y.max()+1
 X=np.concatenate([X,Xn],axis=0)
 Y=np.concatenate([Y,Yn],axis=0)
 return X,Y

И попробуем его кластеризовать:

from sklearn.datasets import make_moons
S,N=1000,5
X,Y=MakeCircles()
fig,axs=plt.subplots(1,3,figsize=(7,2.5),constrained_layout=True)
axs[0].set_title('source clusters')
axs[0].scatter(X[:,0],X[:,1],c=Y,s=1)
arr=[]
for n_cl in range(1,30):
 clf=GaussianMixture(n_components=n_cl)
 clf.fit(X)
 bic=clf.bic(X)
 arr.append([n_cl,bic])
arr=np.array(arr)
axs[1].set_title('BIC')
axs[1].set_xlabel('n_components')
axs[1].set_ylabel('BIC')
axs[1].plot(arr[:,0],arr[:,1])
n_cl=arr[np.argmin(arr[:,1]),0]
clf=GaussianMixture(n_components=int(n_cl))
Yp=clf.fit_predict(X)
axs[2].set_title('BIC-optimal clusters')
axs[2].scatter(X[:,0],X[:,1],c=Yp,s=1)

Получим очень странную кластеризацию:

Рис.6 Случай GM+BIC - вложенные кольца с шумом

Из-за того, что форма кластеров (кольца) отличается от модели (эллипсоиды), у нас получилось больше 20 кластеров.

Не лучше результат и у DBSCAN+Silhouette:

from sklearn.datasets import make_moons
import matplotlib.pyplot as plt
import numpy as np
from sklearn.metrics import silhouette_score
from sklearn.cluster import DBSCAN 
S,N=1000,5
X,Y=MakeCircles()
fig,axs=plt.subplots(1,3,figsize=(7,2.5),constrained_layout=True)
axs[0].set_title('source clusters')
axs[0].scatter(X[:,0],X[:,1],c=Y,s=1)
arr=[]
for eps in np.linspace(0.1,10,100):
 clf=DBSCAN(eps=eps)
 Yp=clf.fit_predict(X)
 if np.unique(Yp).shape[0]>1:
  ssc=silhouette_score(X,Yp)
  arr.append([eps,ssc])
arr=np.array(arr)
axs[1].set_title('Силуэт')
axs[1].set_xlabel('eps')
axs[1].set_ylabel('Silhouette')
axs[1].plot(arr[:,0],arr[:,1])
eps=arr[np.argmax(arr[:,1]),0]
clf=DBSCAN(eps=eps)
Yp=clf.fit_predict(X)
axs[2].set_title('Silh-optimal clusters')
axs[2].scatter(X[:,0],X[:,1],c=Yp,s=1);

Но здесь причина проблемы - в шуме и в плохой работе Силуэта на вложенных кластерах (межкластерное расстояние значительно меньше, чем их собственный размер):

Рис.7 Случай DBSCAN+Silhouette - вложенные кольца с шумом

Напрашивается неутешительный вывод: если ваши кластера сложной формы да еще и в шумах, вам могут не помочь ни DBSCAN+Silhouette, ни GM+BIC. А что делать, если ваши данные именно такие? Оказывается, в этом случае можно попытаться объединить GM и DBSCAN и получить неплохой результат.

Суперкластеризация методом GMsDB

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

Интересной задачей является построение алгоритма кластеризации, который будет с одной стороны устойчив к шуму (как GM), а с другой стороны – сможет разделять кластера сложной формы (как DBSCAN).

Идея такого метода достаточно проста и состоит в последовательном использовании двух этих методов: сначала мы разделим наши данные на кластера с помощью GM, а затем, рассматривая каждый GM кластер, как точку в некотором новом пространстве, объединим их в суперкластера с помощью DBSCAN.

Сложностями в этом случае будут следующие:

как выбрать гиперпараметр у GM;
как выбрать новое пространство, в котором каждый GM кластер будет точкой и как ввести расстояние между такими точками;
как выбрать гиперпараметр у DBSCAN.

Первая задача решается очень просто – использовать критерий BIC. Из рис.6 очевидно, что в случае сложных кластеров мы получим после этого кластеров больше, чем нам нужно, и их позже можно объединить в суперкластера.

Вторая задача несколько сложнее, но решается с помощью известного расстояния Махаланобиса. Расстояние Махаланобиса (https://ru.wikipedia.org/wiki/Расстояние_Махаланобиса) можно рассматривать, как расстояние от точки до центра некого эллипсоида, деленное на ширину эллипсоида в направлении этой точки. Для тех, кто знаком с теорией проверки статистических гипотез это расстояние имеет очень интересную интерпретацию. В одномерном случае, 95% точек эллипсоида расположены на расстоянии Махаланобиса меньшем 1.96. Это означает, что расстояние Махаланобиса можно использовать в качестве статистической меры того, что точка принадлежит заданному кластеру – если расстояние от нее до кластера ниже заданного порога – то точка ему принадлежит, если выше – то нет. На основе этого принципа можно ввести расстояние между кластерами – построить распределение расстояний Махаланобиса между всеми парами точек и взять нижнюю границу 95% доверительного интервала – это и будет расстояние между кластерами. Расстояние это будет неевклидовым, поскольку связано не с декартовым расстоянием между точками, а с вероятностями принадлежности этих точек гауссовым кластерам. Симметризовав это расстояние (взяв из двух расстояний между двумя эллипсоидами наибольшее), можно построить из этих расстояний матрицу расстояний и попробовать объединить исходные гауссовы кластера точек в суперкластера методом DBSCAN. Таким образом, мы перешли в новое пространство, использовав однако не координаты (мы их не знаем), но расстояния между точками. Каждая точка в этом пространстве - это целый GM кластер, а расстояние между ними - зависит обратно вероятности, что два кластера пересекаются (по смыслу расстояния Махаланобиса) - чем меньше вероятность пересечения GM кластеров, тем расстояние больше.

Третья задача еще сложнее. Как мы уже показали – каждая метрика качества зависит от формы кластеров. Поскольку расстояние неевклидово – мы не можем использовать Silhouette, а поскольку исходные кластера гауссовы и построены на основе критерия BIC – то и критерий BIC заново использовать будет не очень разумно. Поэтому нужен новый критерий, основанный на вероятностных методах. Тем более, что и расстояние в нашем пространстве имеет вполне прозрачный смысл вероятности.

Здесь мы можем снова вернуться к магическому числу 1.96, определяющему граничное расстояние, на котором два кластера можно считать соприкасающимися. Это волшебное число зависит от двух величин:

- сколько процентов данных должно быть внутри доверительного интервала (обычно оно определяется статистическим уровенем значимости \alpha, традиционно в статистике принимаемому равным 0.05, а те самые 95% точек = 1 - \alpha);

- какова размерность исходного пространства D (сколько координат используется для описания объекта).

Если учесть эти две величины, то волшебный порог R для расстояний Махаланобиса при заданных \alpha,D можно расчитать по сравнительно несложной формуле через квантили распределения хи-квадрат:

$R=\sqrt{2\cdot Q_{1-\alpha}}$ $Q_{p} : \int_{-\infty}^{Q_p}\chi^2(x,df=D)dx=p$

Где Q - соответствующий (1-\alpha) квантиль распределения хи-квадрат с числом степеней свободы D.

Для объединения кластеров в суперкластера будем использовать DBSCAN, а в качестве матрицы расстояний – матрицу межкластерных расстояний Махалонобиса, посчитанную на предыдущем шаге. Перебирая гиперпараметр DBSCAN (расстояние \eps) от наименьшего (дающего максимальное число суперкластеров) до наибольшего (дающего только один суперкластер) мы будем получать различные количества суперкластеров. На каком значении гиперпараметра DBSCAN остановиться? Лучше всего остановиться тогда, когда мы нашли максимальное число непересекающихся суперкластеров.

Если предположить, что мы каким-то образом объединили кластера в суперкластера, то для того, чтобы они не пересекались, все они должны быть попарно на расстояниях больших R, вычисленного нами для заданного уровня значимости \alpha и размерности данных D. Поскольку для расстояния отдельной точки от эллипса стандартный уровень значимости \alpha=0.05, а мы рассчитываем симметричное расстояние между всеми точками двух эллипсов, имеет смысл выбрать уровень значимости \alpha=2*0.05=0.1 .

Такой критерий был назван матричным критерием (MC):

MC=количество непересекающихся суперкластеров / количество суперкластеров

Если он равен 1, то все суперкластера в нашем разбиении не пересекаются, и мы нашли то, что искали. При этом очевидно, что надо искать максимальное число суперкласетров, при котором выполняется условие MC=1.

Этот алгоритм и реализован в виде класса GMsDB.

Код и результат работы алгоритма показан ниже (требует установки пакета gmsdb через pip).

import matplotlib.pyplot as plt
import numpy as np
from gmsdb import GMSDB 
X,Y=MakeCircles()
fig,axs=plt.subplots(1,2,figsize=(5,2.5),constrained_layout=True)
axs[0].set_title('source clusters')
axs[0].scatter(X[:,0],X[:,1],c=Y,s=1)
arr=[]
clf=GMSDB(n_components=40)
Yp=clf.fit_predict(X)
axs[1].set_title('GMSDB clusters')
axs[1].scatter(X[:,0],X[:,1],c=Yp,s=1);

Параметр n_components=40 в коде - это максимальное число гауссовых компонент при их поиске на первом этапе. Если оно достаточно большое, то от его дальнейшего увеличения зависит только увеличение времени работы.

К недостаткам метода можно отнести его низкое быстродействие и вероятностный характер получаемых разбиений, но необходимые мне задачи он решает (я занимаюсь в том числе и кластеризацией экспериментальных данных).

Какие плюсы? Прогноз кластера для новых точек не требует нового обучения (у метода реализована функция .predict - метод-то модельный и параметры модели мы определили при фитировании). Он может определять не только кластера точек, но и вероятности того, что точка принадлежит каждому из кластеров (.predict_proba у класса тоже работает - опять-же из-за того, что метод модельный и распределения вероятностей мы вычислили) - то есть это - алгоритм мягкой/нечеткой (fuzzy) кластеризации. И да, у алгоритма фактически нет гиперпараметров: все контролирует только уровень статистической значимости \alpha=0.1, выбранный по умолчанию исходя из статистических соображений. Этот уровень значимости годился для всех датасетов, на которых метод проверялся.

Более детально описание метода и примеры его использования можно посмотреть здесь: https://arxiv.org/abs/2309.02623 , а сам код и комментарии к использованию этого метода можно найти на гитхабе: https://github.com/berng/GMSDB

Спасибо Kandinsky за иллюстрацию

Кластеризация: зачем и как объединять Gaussian Mixture и DBSCAN?