A simple mixture model seems to work pretty well for this problem.
In general, to get a point that minimizes the distance to all other points in a dataset, you can just take the mean. In this case, you want to find a point that minimizes the distance from a subset of concentrated points. If you postulate that a point can either come from the concentrated set of points of interest or from a diffuse set of background points, then this gives a mixture model.
I have included some python code below. The concentrated area is modeled by a high-precision normal distribution and the background point are modeled by either a low-precision normal distribution or a uniform distribution over a bounding box on the dataset (there is a line of code that can be commented out to switch between these options). Also, mixture models can be somewhat unstable, so running the EM algorithm a few times with random initial conditions and choosing the run with the highest log-likelihood gives better results.
If you are actually looking at airplanes, then adding some sort of time dependent dynamics will probably improve your ability to infer the home base immensely.
I would also be wary of Rossimo’s formula because it includes some pretty-strong assumptions about crime distributions.
#the dataset
sdata='''41.892694,-87.670898
42.056048,-88.000488
41.941744,-88.000488
42.072361,-88.209229
42.091933,-87.982635
42.149994,-88.133698
42.171371,-88.286133
42.23241,-88.305359
42.196811,-88.099365
42.189689,-88.188629
42.17646,-88.173523
42.180531,-88.209229
42.18168,-88.187943
42.185496,-88.166656
42.170485,-88.150864
42.150634,-88.140564
42.156743,-88.123741
42.118555,-88.105545
42.121356,-88.112755
42.115499,-88.102112
42.119319,-88.112411
42.118046,-88.110695
42.117791,-88.109322
42.182189,-88.182449
42.194145,-88.183823
42.189057,-88.196182
42.186513,-88.200645
42.180917,-88.197899
42.178881,-88.192062
41.881656,-87.6297
41.875521,-87.6297
41.87872,-87.636566
41.872073,-87.62661
41.868239,-87.634506
41.86875,-87.624893
41.883065,-87.62352
41.881021,-87.619743
41.879998,-87.620087
41.8915,-87.633476
41.875163,-87.620773
41.879125,-87.62558
41.862763,-87.608757
41.858672,-87.607899
41.865192,-87.615795
41.87005,-87.62043
42.073061,-87.973022
42.317241,-88.187256
42.272546,-88.088379
42.244086,-87.890625
42.044512,-88.28064
39.754977,-86.154785
39.754977,-89.648437
41.043369,-85.12207
43.050074,-89.406738
43.082179,-87.912598
42.7281,-84.572754
39.974226,-83.056641
38.888093,-77.01416
39.923692,-75.168457
40.794318,-73.959961
40.877439,-73.146973
40.611086,-73.740234
40.627764,-73.234863
41.784881,-71.367187
42.371988,-70.993652
35.224587,-80.793457
36.753465,-76.069336
39.263361,-76.530762
25.737127,-80.222168
26.644083,-81.958008
30.50223,-87.275391
29.436309,-98.525391
30.217839,-97.844238
29.742023,-95.361328
31.500409,-97.163086
32.691688,-96.877441
32.691688,-97.404785
35.095754,-106.655273
33.425138,-112.104492
32.873244,-117.114258
33.973545,-118.256836
33.681497,-117.905273
33.622982,-117.734985
33.741828,-118.092041
33.64585,-117.861328
33.700707,-118.015137
33.801189,-118.251343
33.513132,-117.740479
32.777244,-117.235107
32.707939,-117.158203
32.703317,-117.268066
32.610821,-117.075806
34.419726,-119.701538
37.750358,-122.431641
37.50673,-122.387695
37.174817,-121.904297
37.157307,-122.321777
37.271492,-122.033386
37.435238,-122.217407
37.687794,-122.415161
37.542025,-122.299805
37.609506,-122.398682
37.544203,-122.0224
37.422151,-122.13501
37.395971,-122.080078
45.485651,-122.739258
47.719463,-122.255859
47.303913,-122.607422
45.176713,-122.167969
39.566,-104.985352
39.124201,-94.614258
35.454518,-97.426758
38.473482,-90.175781
45.021612,-93.251953
42.417881,-83.056641
41.371141,-81.782227
33.791132,-84.331055
30.252543,-90.439453
37.421248,-122.174835
37.47794,-122.181702
37.510628,-122.254486
37.56943,-122.346497
37.593373,-122.384949
37.620571,-122.489319
36.984249,-122.03064
36.553017,-121.893311
36.654442,-121.772461
36.482381,-121.876831
36.15042,-121.651611
36.274518,-121.838379
37.817717,-119.569702
39.31657,-120.140991
38.933041,-119.992676
39.13785,-119.778442
39.108019,-120.239868
38.586082,-121.503296
38.723354,-121.289062
37.878444,-119.437866
37.782994,-119.470825
37.973771,-119.685059
39.001377,-120.17395
40.709076,-73.948975
40.846346,-73.861084
40.780452,-73.959961
40.778829,-73.958931
40.78372,-73.966012
40.783688,-73.965325
40.783692,-73.965615
40.783675,-73.965741
40.783835,-73.965873
'''
import StringIO
import numpy as np
import re
import matplotlib.pyplot as plt
def lp(l):
return map(lambda m: float(m.group()),re.finditer('[^, n]+',l))
data=np.array(map(lp,StringIO.StringIO(sdata)))
xmn=np.min(data[:,0])
xmx=np.max(data[:,0])
ymn=np.min(data[:,1])
ymx=np.max(data[:,1])
# area of the point set bounding box
area=(xmx-xmn)*(ymx-ymn)
M_ITER=100 #maximum number of iterations
THRESH=1e-10 # stopping threshold
def em(x):
print 'nSTART EM'
mlst=[]
mu0=np.mean( data , 0 ) # the sample mean of the data - use this as the mean of the low-precision gaussian
# the mean of the high-precision Gaussian - this is what we are looking for
mu=np.random.rand( 2 )*np.array([xmx-xmn,ymx-ymn])+np.array([xmn,ymn])
lam_lo=.001 # precision of the low-precision Gaussian
lam_hi=.1 # precision of the high-precision Gaussian
prz=np.random.rand( 1 ) # probability of choosing the high-precision Gaussian mixture component
for i in xrange(M_ITER):
mlst.append(mu[:])
l_hi=np.log(prz)+np.log(lam_hi)-.5*lam_hi*np.sum((x-mu)**2,1)
#low-precision normal background distribution
l_lo=np.log(1.0-prz)+np.log(lam_lo)-.5*lam_lo*np.sum((x-mu0)**2,1)
#uncomment for the uniform background distribution
#l_lo=np.log(1.0-prz)-np.log(area)
#expectation step
zs=1.0/(1.0+np.exp(l_lo-l_hi))
#compute bound on the likelihood
lh=np.sum(zs*l_hi+(1.0-zs)*l_lo)
print i,lh
#maximization step
mu=np.sum(zs[:,None]*x,0)/np.sum(zs) #mean
lam_hi=np.sum(zs)/np.sum(zs*.5*np.sum((x-mu)**2,1)) #precision
prz=1.0/(1.0+np.sum(1.0-zs)/np.sum(zs)) #mixure component probability
try:
if np.abs((lh-old_lh)/lh)<THRESH:
break
except:
pass
old_lh=lh
mlst.append(mu[:])
return lh,lam_hi,mlst
if __name__=='__main__':
#repeat the EM algorithm a number of times and get the run with the best log likelihood
mx_prm=em(data)
for i in xrange(4):
prm=em(data)
if prm[0]>mx_prm[0]:
mx_prm=prm
print prm[0]
print mx_prm[0]
lh,lam_hi,mlst=mx_prm
mu=mlst[-1]
print 'best loglikelihood:', lh
#print 'final precision value:', lam_hi
print 'point of interest:', mu
plt.plot(data[:,0],data[:,1],'.b')
for m in mlst:
plt.plot(m[0],m[1],'xr')
plt.show()
Метод
-средних
на примере
Изображение
Описание
Шаг 1.
Определим
число кластеров, на которое требуется
разбить исходное множество
.
Шаг 2.
Случайным образом выберем
две точки, которые будут начальными
центрами кластеров. Пусть это будут
точки
и
.
На рисунке они закрашены черным
цветом.
# |
X |
Y |
Расстояние от |
Расстояние от |
Номер кластера |
A |
1 |
3 |
2.00 |
2.24 |
1 |
B |
3 |
3 |
2.83 |
2.24 |
2 |
C |
4 |
3 |
3.61 |
2.83 |
2 |
D |
5 |
3 |
4.47 |
3.61 |
2 |
E |
1 |
2 |
1.00 |
1.41 |
1 |
F |
4 |
2 |
3.16 |
2.24 |
2 |
G |
1 |
1 |
0.00 |
1.00 |
1 |
H |
2 |
1 |
1.00 |
0.00 |
2 |
Цветом выделено минимальное |
Шаг
3.
Проход
1.
Для
каждой точки определим к ней ближайший
центр кластера с помощью расстояния
Евклида. В таблице представлены
вычисленные с помощью формулы
расстояния
между центрами кластеров
,
и каждой точкой исходного множества,
а также указано, к какому кластеру
принадлежит та или иная точка.
Таким
образом, кластер 1 содержит точки A,
E, G, а кластер 2 – точки B,
C, D, F, H. Как только
определятся члены кластеров, может
быть рассчитана сумма квадратичных
ошибок (
— точки
-того
кластера):
Шаг
4.
Проход
1.
Для
каждого кластера вычисляется его
центроид, и центр кластера перемещается
в него.
Центроид
для кластера 1:
.
Центроид
для кластера 2:
.
Начальные
центры кластеров закрашены зеленым
цветом, а центроиды, вычисленные при
1-м проходе алгоритма, – закрашены
черным цветом. Они и будут являться
новыми центрами кластеров, к которым
будет определяться принадлежность
точек данных на втором проходе.
# |
X |
Y |
Расстояние от |
Расстояние от |
Номер кластера |
A |
1 |
3 |
1.00 |
2.67 |
1 |
B |
3 |
3 |
2.24 |
0.75 |
2 |
C |
4 |
3 |
3.16 |
0.72 |
2 |
D |
5 |
3 |
4.12 |
1.52 |
2 |
E |
1 |
2 |
0.00 |
2.63 |
1 |
F |
4 |
2 |
3.00 |
0.57 |
2 |
G |
1 |
1 |
1.00 |
2.95 |
1 |
H |
2 |
1 |
1.41 |
2.13 |
1 |
Цветом выделено минимальное |
Шаг
3.
Проход
2.
После
того, как найдены новые центры
кластеров, для каждой точки снова
определяется ближайший к ней центр
и ее отношение к соответствующему
кластеру. Для это еще раз вычисляются
евклидовы расстояния между точками
и центрами кластеров.
Относительно
большое изменение
привело к тому, что запись H
оказалась ближе к центру
,
что автоматически сделало ее членом
кластера 1. Все остальные записи
остались в тех же кластерах, что и на
предыдущем проходе алгоритма. Таким
образом, кластер 1 будет A, E,
G, H, а кластер 2 – B, C,
D, F. Новая сумма квадратов
ошибок составит (
— точки
-того
кластера):
Результат
показывает уменьшение ошибки
относительно начального состояния
центров кластеров (которая на первом
проходе составляла 36). Это говорит об
улучшении качества кластеризации,
т.е. более высокую «кучность» объектов
относительно центра кластера.
Шаг
4.
Проход
2.
Для
каждого кластера вновь вычисляется
его центроид, и центр кластера
перемещается в него.
Новый
центроид для 1-го кластера:
.
Центроид
для кластера 2:
.
Расположение
кластеров и центроидов после второго
прохода алгоритма представлено на
рисунке.
По
сравнению с предыдущим проходом
центры кластеров изменились
незначительно.
# |
X |
Y |
Расстояние от |
Расстояние от |
Номер кластера |
A |
1 |
3 |
1.27 |
3.01 |
1 |
B |
3 |
3 |
2.15 |
1.03 |
2 |
C |
4 |
3 |
3.02 |
0.25 |
2 |
D |
5 |
3 |
3.95 |
1.03 |
2 |
E |
1 |
2 |
0.35 |
3.09 |
1 |
F |
4 |
2 |
2.76 |
0.75 |
2 |
G |
1 |
1 |
0.79 |
3.47 |
1 |
H |
2 |
1 |
1.06 |
2.66 |
1 |
Цветом выделено минимальное |
Шаг
3.
Проход
3.
Для
каждой записи вновь ищется ближайший
к ней центр кластера.
Следует
отметить, что записей, сменивших
кластер на данном проходе алгоритма,
не было. Новая сумма квадратов ошибок
составит:
Таким
образом, изменение суммы квадратов
ошибок является незначительным по
сравнению с предыдущим проходом.
Шаг
4.
Проход
3.
Для
каждого кластера вновь вычисляется
его центроид, и центр кластера
перемещается в него. Но поскольку на
данном проходе ни одна запись не
изменила своего членства в кластерах,
то положение центроида не поменялось,
и алгоритм завершает свою работу.
Блок-схема
Задача кластеризации
В задаче классификации мы имели дело с восстановлением отображения из множества объектов в конечный набор меток классов. При этом классы были зафиксированы заранее, то есть мы с самого начала примерно понимали, какого рода объекты должны относиться к каждому из них, и мы располагали обучающей выборкой с примерами объектов и классов, к которым они относятся. В задаче кластеризации мы тоже разбиваем объекты на конечное множество классов, но у нас нет ни обучающей выборки, ни понимания, какой будет природа этих классов. То, что модель кластеризации какие-то объекты сочла «похожими», отнеся к одному классу, будет новой информацией, «открытием», сделанным этой моделью. Обучающей выборки у нас также не будет: ведь мы не знаем заранее, что за классы получатся (а иногда и сколько их будет). Таким образом, кластеризация — это задача обучения без учителя. Из-за общего сходства постановок задач в литературе кластеризацию иногда называют unsupervised classification.
Методы кластеризации часто применяют, когда фактически нужно решить задачу классификации, но обучающую выборку собрать затруднительно (дорого или долго). При этом валидационную выборку для оценки результатов кластеризации собрать значительно проще, так как для неё требуется меньше примеров. При этом стоит помнить, что точность работы supervised-методов значительно выше. Поэтому, если обучающую выборку всё-таки можно собрать, лучше решать задачу классификации, чем задачу кластеризации.
Примеры задач кластеризации
Хороший пример применения методов кластеризации — анализ геоданных. В мобильных приложениях, собирающих геоданные пользователей, часто требуется понять, где именно пользователь находился. GPS-координаты известны с некоторой погрешностью, пользователь тоже обычно двигается, поэтому вместо точного положения часто приходится иметь дело с роем точек. Положение усугубляется, когда мы пытаемся анализировать поведение сразу тысяч людей в какой-то локации — например, определить, в каких точках люди чаще всего садятся в такси у аэропорта. Может показаться, что достаточно посмотреть на данные — и мы увидим в точности нужные нам кластеры. Изображение ниже показывает, как может выглядеть ситуация всего для нескольких пользователей: согласно данным GPS, такси подбирают пассажиров и внутри здания аэропорта, и на взлётной полосе, и там, где это происходит на самом деле:
Подобная задача решалась в Яндекс.Такси при разработке пикап-пойнтов (наиболее удобных точек вызова такси, подсвечиваемых в приложении). Координаты точек заказа кластеризовались таким образом, чтобы кластер соответствовал какому-то одному, удобному для пользователя месту, и центры кластеров использовались как кандидаты в пикап-пойнты. Те кандидаты, которые удовлетворяли простым фильтрам (например, не попадали в здание или в воду), использовались в приложении. При этом не обходилось и без вручную проставленных пикап-пойнтов: например, такое решение использовалось в окрестностях аэропортов.
Другой пример кластеризации геоданных, который всегда рядом с нами, — это интерфейсы для просмотра фотографий в вашем смартфоне. Почти наверняка вы можете просмотреть их в привязке к местам, где они были сделаны, и по мере масштабирования карты вы будете видеть разное количество кластеров фотографий. Кстати, если говорить об интерфейсах, то есть и другой интересный пример: если нужно подстроить цветовую схему вашего интерфейса под выбираемое пользователем изображение (например, фоновую картинку), достаточно кластеризовать цвета из пользовательского изображения, используя RGB-представление (или любое другое) как признаки цвета, и воспользоваться для оформления цветами, соответствующими центрам кластеров.
Простейшие методы кластеризации с помощью графов
Можно приводить примеры не только про геоаналитику, однако тема геоданных поможет нам придумать пару наиболее простых и наглядных методов кластеризации. Представим, что перед нами рой геокоординат и нам нужно предложить по этим данным пикап-пойнты для такси. Разберём пару очевидных методов.
Выделение компонент связности
Логично попробовать объединить точки, которые находятся друг от друга на расстоянии двух-трёх метров, а потом просто выбрать наиболее популярные места. Для этого давайте построим на известных нам точках граф: точки, расстояние между которыми в пределах трёх метров, мы соединим рёбрами. Выделим в этом графе компоненты связности, они и будут нашими кластерами.
У этого способа есть пара очевидных недостатков. Во-первых, может найтись сколько угодно длинная цепочка точек, в которой соседние отстоят друг от друга на пару метров, — и вся она попадёт в одну компоненту связности. В итоге наша отсечка по трём метрам имеет очень опосредованное отношение к диаметру кластеров, а сами кластеры будут получаться значительно больше, чем нам хотелось бы. Во-вторых (и с первой проблемой это тоже связано), непонятно, как мы выбираем максимальное расстояние, при котором соединяем точки ребром. В данной задаче ещё можно предъявить хоть какую-то логику, а вот если бы мы кластеризовали не геометки, а что-то многомерное, например электронные письма по их тематике, придумать отсечку было бы уже сложнее. Если наша цель — не только решить практическую задачу, но и придумать достаточно общий метод кластеризации, понятно, что нам хочется понимать, как подбирать параметры этого метода (в данном случае условие добавления рёбер в граф). Эти соображения могут привести нас к другому решению.
Минимальное остовное дерево
Вместо того чтобы проводить рёбра в графе, давайте их удалять. Построим минимальное остовное дерево, считая расстояния между точками весами рёбер. Тогда, удалив $N$ рёбер с наибольшим весом, мы получим $N+1$ компоненту связности, которые, как и в прошлом подходе, будем считать кластерами. Различие в том, что теперь нам нужно задавать не расстояние, при котором проводится ребро, а количество кластеров. С одной стороны, если мы решаем задачу расчёта пикап-пойнтов в какой-то конкретной локации (аэропорт, торговый центр, жилой дом), нам может быть понятно, сколько примерно пикап-пойнтов мы хотим получить. С другой стороны, даже без локального рассмотрения можно просто сделать достаточно много кластеров, чтобы было из чего выбирать, но при этом достаточно мало, чтобы в каждый кластер попадало репрезентативное количество точек. Аналогичная логика будет справедлива и во многих других задачах кластеризации: количество кластеров — достаточно общий и достаточно хорошо интерпретируемый параметр, чтобы настраивать его вручную, поэтому во многих методах кластеризации количество кластеров выступает как гиперпараметр.
Далее будем рассматривать некоторую обобщённую задачу кластеризации без привязки к нашему примеру с анализом геоданных. Мы приведём три наиболее популярных метода кластеризации — k-средних, иерархическую кластеризацию и DBSCAN, а затем рассмотрим вопросы оценки качества кластеризации.
Метод K средних
Пожалуй, один из наиболее популярных методов кластеризации — это метод K-средних (K-means). Основная идея метода — итеративное повторение двух шагов:
- распределение объектов выборки по кластерам;
- пересчёт центров кластеров.
В начале работы алгоритма выбираются $K$ случайных центров в пространстве признаков. Каждый объект выборки относят к тому кластеру, к центру которого объект оказался ближе. Далее центры кластеров пересчитывают как среднее арифметическое векторов признаков всех вошедших в этот кластер объектов (то есть центр масс кластера). Как только мы обновили центры кластеров, объекты заново перераспределяются по ним, а затем можно снова уточнить положение центров. Процесс продолжается до тех пор, пока центры кластеров не перестанут меняться.
Выбор начального приближения
Первый вопрос при выборе начального положения центров — как, выбирая центры из некоторого случайного распределения, не попасть в область пространства признаков, где нет точек выборки. Базовое решение — просто выбрать в качестве центров какие-то из объектов выборки.
Вторая потенциальная проблема — кучное размещение центров. В этом случае их начальное положение с большой вероятностью окажется далёким от итогового положения центров кластеров. Например, для таких изначальных положений центров
мы получим неправильную кластеризацию.
Чтобы бороться с этим явлением, выгодно брать максимально удаленные друг от друга центры.
На практике работает следующая эвристика:
-
первый центр выбираем случайно из равномерного распределения на точках выборки;
-
каждый следующий центр выбираем из случайного распределения на объектах выборки, в котором вероятность выбрать объект пропорциональна квадрату расстояния от него до ближайшего к нему центра кластера. Модификация K-means, использующая эту эвристику для выбора начальных приближений, называется K-means++.
Выбор метрик
Так как работа метода K-средних состоит из последовательного повторения до сходимости двух шагов, обоснованность применения различных метрик (расстояний между точками, а не метрик качества 🙂 или функций близости связана с тем, «ломают» они какой-либо из этих шагов или нет.
Первый шаг с отнесением объектов к ближайшим центрам не зависит от вида метрики. Второй шаг предполагает пересчёт центров как среднего арифметического входящих в кластер точек, и вот здесь будет подвох: к оптимальности выбора центров в среднем арифметическом приводит именно евклидова метрика (подробнее в разделе «Что оптимизирует K-means»).
Однако на практике никто не мешает использовать метод и без должного обоснования, поэтому можно экспериментировать с любыми расстояниями, с той лишь оговоркой, что не будет никаких теоретических гарантий, что метод сработает. Наиболее распространённая альтернатива евклидовой метрике — это косинусная мера близости векторов (она особенно популярна в задачах анализа текстов):
$$CosineSimilarity(mu_k, x_i)=frac{<mu_k, x_i>}{|mu_k|_2 cdot |x_i|_2}$$
При её использовании стоит не забывать, что косинусная мера — это функция близости, а не расстояние, так что чем больше её значения, тем ближе друг к другу векторы.
Mini-batch K-means
Несложно заметить, что, если считать $K$ и размерность пространства признаков константами, оба шага алгоритма работают за $O(n)$, где n — количество объектов обучающей выборки. Отсюда возникает идея ускорения работы алгоритма. В mini-batch K-means мы не считаем шаги сразу на всей выборке, а на каждой итерации выбираем случайную подвыборку (мини-батч) и работаем на ней. В случае когда исходная выборка очень велика, переход к пакетной обработке не приводит к большой потере качества, зато значительно ускоряет работу алгоритма.
Понижение размерности
С другой стороны, вычисление расстояний и средних делается за $O(d)$, где $d$ — размерность пространства признаков, так что другая идея ускорения K-means — это предварительно понизить размерность пространства признаков (с помощью PCA или эмбеддингов). Особенно удачно эта идея работает в задачах кластеризации текстов, когда K-means применяют на эмбеддингах слов: получается выиграть не только в скорости работы, но и в интерпретируемости результатов кластеризации.
Кстати, сам алгоритм кластеризации тоже можно использовать как метод понижения размерности. Если вы решаете задачу обучения с учителем и пространство признаков очень разнообразно (то есть обучающая выборка не даёт вам достаточно статистики при столь большом числе признаков), можно выполнить кластеризацию объектов выборки на 500 или 1000 кластеров и оперировать попаданием объектов в какой-то кластер как признаком. Такой подход называется квантизацией пространства признаков (feature space quantization) и часто помогает на практике, когда нужно огрубить признаки, добавить им интерпретируемости или же, наоборот, обезличить.
Хрестоматийный пример такого использования кластеризации — метод bag of visual words, расширяющий bag of words из анализа текстов на работу с изображениями. Идея метода в том, чтобы строить признаковое описание изображений на основе входящих в него фрагментов: так, изображения с лицами будут содержать фрагменты с носом, глазами, ртом, а изображения с машинами — колёса, зеркала, двери. Но проблема здесь в том, что нарезать такие фрагменты из обучающей выборки и искать точные совпадения в новых примерах изображений, которые нужно классифицировать, — безнадёжная затея. В жизни фрагменты изображений не повторяются в других изображениях с попиксельной точностью. Решение этой проблемы оказалось возможным при помощи алгоритмов кластеризации (исторически использовался именно K-means): фрагменты изображений из обучающей выборки кластеризовали на 100–1000 кластеров («визуальных слов»), а проходясь по новым изображениям, также нарезали их на фрагменты и относили к одному из этих кластеров. В итоге как новые изображения, так и изображения из обучающей выборки можно было описать количеством вхождений в них фрагментов из различных кластеров («визуальных слов»), так же как в анализе текстов описывают текст количеством вхождений в него слов из словаря. В таком признаковом пространстве уже можно было успешно обучать модели машинного обучения.
Сейчас выделение «визуальных слов» в задаче классификации изображений происходит автоматически: с одной стороны, задачи компьютерного зрения теперь решаются нейросетями, но с другой стороны — если мы визуализируем отдельные слои этих нейросетей, станет понятно, что их логика работы во многом похожа на описанную выше. При этом идея квантизации признаков не утратила своей актуальности. Вот лишь несколько современных примеров её применения:
- Если вам необходимо дать возможность заказчику (например, внешней компании) анализировать используемые вами признаки — отсутствие провалов в данных и какие-то другие общие показатели, но нельзя отдавать признаки как есть (например, из-за законодательства, регулирующего передачу пользовательских данных), возможное решение — это агрегировать признаки по кластерам.
- Та же цель может быть отчасти достигнута, если перейти к самим кластерам как к признакам, чтобы скрыть исходные признаки.
- Переход к кластерам может быть сделан не с целью что-то скрыть, а наоборот, с целью повысить интерпретируемость: исходные сырые данные часто не вполне понятны бизнесу, но позволяют построить маркетинговые сегменты по различным коммерческим интересам пользователей, из-за чего становится удобно показывать принадлежность к этим сегментам как исходные признаки, не вдаваясь в детали о том, на каких данных эти сегменты построены.
- Для ускорения поиска похожих объектов в пространстве признаков вы также можете проводить поиск внутри того же кластера и соседних кластеров, так что за счёт «огрублённого» вида признаков какие-то процессы можно ещё и ускорить.
Что оптимизирует K-means
Проговорим на интуитивном уровне, какую оптимизационную задачу решает K-means.
Оба шага алгоритма работают на уменьшение среднего квадрата евклидова расстояния от объектов до центров их кластеров:
$$Phi_0 = frac{1}{nK} sumlimits_{k=1}^{K} sumlimits_{i=1}^{n} (mu_k – x_i)^2 mathbb{I}[a(x_i)=k]$$
На шаге отнесения объектов к одному из кластеров мы выбираем кластер с ближайшим центроидом, то есть минимизируем каждое слагаемое в $Phi_0$: все потенциально большие слагаемые мы зануляем, а оставляем ненулевыми только наименьшие из возможных (при условии фиксирования центров кластеров).
На шаге пересчёта центров кластеров мы выбираем центр таким образом, чтобы при фиксированном наборе объектов, относящихся к кластеру, для всех $k$ минимизировать выражение, стоящее под суммой по $k$:
$$sumlimits_{i=1}^{n} (mu_k – x_i)^2 mathbb{I}[a(x_i)=k]$$
Здесь уже становится принципиально, что мы определяем квадрат расстояния как квадрат разности векторов, так как именно отсюда при дифференцировании по $mu_k$ и записи необходимого условия экстремума получается, что центры кластеров нужно пересчитывать как средние арифметические $x_i$, принадлежащих кластеру.
Этих соображений, конечно, недостаточно, чтобы утверждать, что мы найдём минимум $Phi_0$. Более того, гарантии того, что мы найдём глобальный минимум, вообще говоря, нет. Однако, потратив чуть больше усилий, можно доказать, что процесс сойдётся в один из локальных минимумов.
Также можно справедливо заметить, что, так как любой центр кластера — это среднее арифметическое входящих в кластер объектов $x_i$, на выборке фиксированного размера есть лишь конечное множество потенциальных центров кластеров. Если предположить, что в ходе работы K-means не зацикливается, отсюда следует, что рано или поздно центры кластеров не изменятся на следующем шаге и алгоритм сойдётся. При этом фактическая сходимость, конечно же, происходит задолго до полного перебора всех возможных центров кластеров.
Иерархическая агломеративная кластеризация
Другой классический метод кластеризации — это иерархическая кластеризация. Иногда дополнительно уточняют: иерархическая агломеративная кластеризация. Название указывает сразу на два обстоятельства.
Во-первых, есть деление алгоритмов кластеризации на агломеративные (agglomerative) и дивизивные, или дивизионные (divisive). Агломеративные алгоритмы начинают с небольших кластеров (обычно с кластеров, состоящих из одного объекта) и постепенно объединяют их в кластеры побольше. Дивизивные начинают с больших кластеров (обычно – с одного единственного кластера) и постепенно делят на кластеры поменьше.
Во-вторых, кластеризация бывает, по аналогии с оргструктурой в организациях, плоской (когда все кластеры равноправны и находятся на одном уровне кластеризации) и иерархической (когда кластеры бывают вложены друг в друга и образуют древовидную структуру).
В случае иерархической агломеративной кластеризации мы действительно будем начинать с кластеров из одного объекта, постепенно объединяя их, а уже последовательность этих объединений даст структуру вложенности кластеров. Даже если в итоге мы будем использовать кластеры с одного уровня, не углубляясь ни в какую вложенность, кластеризация всё равно называется иерархической, так как иерархия естественным образом возникает в процессе работы алгоритма.
Сам алгоритм выглядит предельно просто:
- Создаём столько кластеров, сколько у нас объектов в выборке, каждый объект — в своём отдельном кластере.
- Повторяем итеративно слияние двух ближайших кластеров, пока не выполнится критерий останова.
Расстояния в иерархической кластеризации
Как измерить расстояние между кластерами из одного объекта? Нужно просто взять расстояние между этими объектами. Остаётся вопрос, как обобщить расстояние между объектами до расстояния между кластерами (если в них более одного объекта). Традиционные решения — брать среднее расстояние между объектами кластеров, минимальное расстояние или максимальное. Если обозначить кластеры $U$ и $V$, расстояние между ними в этом случае можем вычислять по одной из формул:
$$ d_{avg}(U, V) = frac{1}{|U| cdot |V|} sumlimits_{u in U} sumlimits_{v in V} rho(u,v) $$
$$ d_{min}(U, V) = minlimits_{(u,v) in U times V} rho(u,v) $$
$$ d_{max}(U, V) = maxlimits_{(u,v) in U times V} rho(u,v) $$
Используемая формула расстояния между кластерами — один из гиперпараметров алгоритма. Кроме приведённых стандартных вариантов бывают и более экзотичные, например расстояние Уорда (Ward distance). В наиболее общем виде способы задания расстояния между кластерами даются формулой Ланса — Уильямса (Lance — Williams; более подробно вы можете почитать в этой статье). Сами же расстояния между объектами можно задавать любой метрикой, как евклидовой, так и манхэттенским расстоянием или, например, косинусной мерой (с той лишь поправкой, что это мера близости, а не расстояние).
Условия окончания работы алгоритма (критерии останова)
В качестве условия для завершения работы алгоритма можем выбрать либо получение нужного количества кластеров (количество кластеров может быть гиперпараметром алгоритма), либо выполнение эвристик на основе расстояния между объединяемыми кластерами (например, если расстояние сливаемых кластеров значительно выросло по сравнению с прошлой итерацией). На практике же обычно кластеризацию проводят вплоть до одного кластера, включающего в себя всю выборку, а затем анализируют получившуюся иерархическую структуру с помощью дендрограммы.
Дендрограмма
По мере объединения кластеров, каждой итерации алгоритма соответствует пара объединяемых на этой итерации кластеров, а также расстояние между кластерами в момент слияния. Расстояния с ростом итерации будут только увеличиваться, поэтому возникает возможность построить следующую схему, называемую дендрограммой:
Здесь по горизонтали внизу отмечены объекты кластеризуемой выборки, под горизонтальной осью подписаны номера объектов, а их расположение вдоль оси продиктовано только эстетическими соображениями: нам удобно строить дендрограмму так, чтобы никакие дуги в ней не пересекались. По вертикали отложены расстояния между кластерами в момент слияния. Когда происходит объединение кластеров, состоящих из нескольких объектов, соответствующая этой итерации алгоритма дуга идёт не до конкретных объектов выборки, а до дуги другого кластера.
Таким образом мы получаем наглядную визуализацию древовидной структуры процесса кластеризации. В частности, на дендрограмме мы можем визуально заметить, в какой момент происходит скачок расстояний между кластерами, и попытаться определить «естественное» количество кластеров в нашей задаче. На практике же это соображение, как правило, остаётся лишь красивой теорией, так как любую кластеризацию можно делать в разной степени «мелкой» и «естественного» количества кластеров в практических задачах часто не существует. В случае же если данные были получены таким образом, что в них действительно есть какое-то естественное количество кластеров, иерархическая кластеризация обычно справляется с определением числа кластеров по дендрограмме заметно хуже, чем DBSCAN. Именно алгоритму DBSCAN мы и посвятим следующий раздел.
DBSCAN
Алгоритм DBSCAN (Density-based spatial clustering of applications with noise) развивает идею кластеризации с помощью выделения связных компонент.
Прежде чем перейти к построению графа, введём понятие плотности объектов выборки в пространстве признаков. Плотность в DBSCAN определяется в окрестности каждого объекта выборки $x_i$ как количество других точек выборки в шаре $B(varepsilon, x_i)$. Кроме радиуса $varepsilon$ окрестности в качестве гиперпараметра алгоритма задается порог $N_0$ по количеству точек в окрестности.
Далее все объекты выборки делятся на три типа: внутренние / основные точки (core points), граничные (border points) и шумовые точки (noise points). К основным относятся точки, в окрестности которых больше $N_0$ объектов выборки. К граничным — точки, в окрестности которых есть основные, но общее количество точек в окрестности меньше $N_0$. Шумовыми называют точки, в окрестности которых нет основных точек и в целом содержится менее $N_0$ объектов выборки.
Алгоритм кластеризации выглядит следующим образом:
- Шумовые точки убираются из рассмотрения и не приписываются ни к какому кластеру.
- Основные точки, у которых есть общая окрестность, соединяются ребром.
- В полученном графе выделяются компоненты связности.
- Каждая граничная точка относится к тому кластеру, в который попала ближайшая к ней основная точка.
Удобство DBSCAN заключается в том, что он сам определяет количество кластеров (по модулю задания других гиперпараметров — $varepsilon$ и $N_0$), а также в том, что метод успешно справляется даже с достаточно сложными формами кластеров. Кластеры могут иметь вид протяжённых лент или быть вложенными друг в друга как концентрические гиперсферы. На изображении ниже показан пример выделения кластеров достаточно сложной формы с помощью DBSCAN:
DBSCAN — один из самых сильных алгоритмов кластеризации, но работает он, как правило, заметно дольше, чем mini-batch K-means, к тому же весьма чувствителен к размерности пространства признаков, поэтому используется на практике DBSCAN только тогда, когда успевает отрабатывать за приемлемое время.
Какой метод кластеризации выбирать?
Если сравнивать частоту использования K-means, иерархической кластеризации и DBSCAN, то на первом месте, бесспорно, будет K-means, а второе место будут делить иерархический подход и DBSCAN. Иерархическая кластеризация — более известный и простой в понимании метод, чем DBSCAN, но довольно редко отрабатывающий качественно. Частая проблема иерархической кластеризации — раннее образование одного гигантского кластера и ряда очень небольших, что приводит к сильной несбалансированности количества объектов в итоговых кластерах. В то же время DBSCAN — менее широко известный подход, но, когда его можно применить, качество, как правило, получается выше, чем в K-means или иерархической кластеризации.
Оценка качества кластеризации
Далее приведём список основных метрик качества кластеризации и обсудим некоторые особенности их применения.
Среднее внутрикластерное растояние
Смысл среднего внутрикластерного расстояния максимально соответствует названию:
$$F_0 = frac{sumlimits_{i=1}^{n} sumlimits_{j=i}^{n} rho(x_i, x_j) mathbb{I}[a(x_i)=a(x_j)]}{sumlimits_{i=1}^{n} sumlimits_{j=i}^{n} mathbb{I}[a(x_i)=a(x_j)]} $$
Сумма расстояний между точками из одного и того же кластера делится на количество пар точек, принадлежащих к одному кластеру. В приведённой выше формуле пары вида $(x_i, x_i)$ включены в рассмотрение, чтобы избежать неопределённости $frac{0}{0}$ в случае, когда в каждом кластере ровно по одному объекту. Однако иногда записывают суммы по $i < j$, просто доопределяя $F_0$ в описанном случае нулём.
Решая задачу кластеризации, мы хотим по возможности получать как можно более кучные кластеры, то есть минимизировать $F_0$.
В случае если у кластеров есть центры $mu_k$, часто рассматривается аналогичная метрика — средний квадрат внутрикластерного расстояния:
$$Phi_0 = frac{1}{nK} sumlimits_{k=1}^{K} sumlimits_{i=1}^{n} rho(mu_k, x_i)^2 mathbb{I}[a(x_i)=k]$$
Среднее межкластерное расстояние
Аналогично среднему внутрикластерному расстоянию вводится среднее межкластерное:
$$F_1 = frac{sumlimits_{i=1}^{n} sumlimits_{j=i}^{n} rho(x_i, x_j) mathbb{I}[a(x_i) neq a(x_j)]}{sumlimits_{i=1}^{n} sumlimits_{j=i}^{n} mathbb{I}[a(x_i) neq a(x_j)]} $$
Среднее межкластерное расстояние, напротив, нужно максимизировать, то есть имеет смысл выделять в разные кластеры наиболее удалённые друг от друга объекты.
Гомогенность
Для измерения следующих метрик (гомогенности, полноты и V-меры) нам уже потребуется разметка выборки. Будем обозначать кластеры, к которым наш алгоритм кластеризации относит каждый объект, буквами $k in lbrace 1, …, K rbrace$, а классы, к которым объекты отнесены разметкой, — буквами $с in lbrace 1, …, С rbrace$. Разумный вопрос при наличии разметки — зачем нам решать задачу кластеризации, если с разметкой можно поставить задачу как задачу классификации. Это и правда хороший вопрос в том случае, если размеченных данных достаточно много для обучения классификатора. На практике же часто встречаются ситуации, когда данных достаточно для оценки качества кластеризации, но всё ещё не хватает для использования методов обучения с учителем.
Пусть $n$ — общее количество объектов в выборке, $n_k$ — количество объектов в кластере номер $k$, $m_c$ — количество объектов в классе номер $с$, а $n_{ck}$ — количество объектов из класса $c$ в кластере $k$. Рассмотрим следующие величины:
$$ H_{class} = – sumlimits_{c = 1}^{C} frac{m_c}{n} log {frac{m_c}{n}} $$
$$ H_{clust} = – sumlimits_{k = 1}^{K} frac{n_k}{n} log {frac{n_k}{n}} $$
$$ H_{class | clust} = – sumlimits_{c = 1}^{C}sumlimits_{k = 1}^{K} frac{n_{ck}}{n_k} log {frac{n_{ck}}{n_k}} $$
Несложно заметить, что эти величины соответствуют формуле энтропии и условной энтропии для мультиномиальных распределений $frac{m_c}{n}$, $frac{n_k}{n}$ и $frac{n_{ck}}{n_k}$ соответственно.
Гомогенность кластеризации определяется следующим выражением:
$$ Homogeneity = 1 – frac{H_{class | clust}}{H_{class}} $$
Отношение $ frac{H_{class vert clust}}{H_{class}} $ показывает, во сколько раз энтропия изменяется за счёт того, что мы считаем известной принадлежность объектов к выделенным нашим алгоритмом кластерам. Худший случай — когда отношение оказывается равным единице (энтропия не изменилась, условная энтропия совпала с обычной), лучший — когда каждый кластер содержит элементы только одного класса и номер кластера, таким образом, точно определяет номер класса (в этом случае $h=1$).
Тривиальный способ максимизировать гомогенность кластеризации — выделить каждый объект выборки в отдельный кластер.
Полнота
Полнота задаётся аналогично гомогенности, с той лишь разницей, что вводится величина $ H_{clust vert class} $, симметричная $ H_{class vert clust} $:
$$ Completeness = 1 – frac{H_{clust | class}}{H_{clust}} $$
Полнота равна единице, когда все объекты класса всегда оказываются в одном кластере.
Тривиальный способ максимизировать полноту кластеризации — объединить всю выборку в один кластер.
V-мера
Гомогенность и полнота кластеризации – это в некотором смысле аналоги точности и полноты классификации. Аналог F-меры для задачи кластеризации тоже есть, он называется V-мерой и связан с гомогенностью и полнотой той же формулой, что и F-мера с точностью и полнотой:
$$ V_{beta} = frac{(1 + beta) cdot Homogeneity cdot Completeness}{beta cdot Homogeneity + Completeness}$$
В частности, $V_1$ по аналогии с $F_1$-мерой в классификации (не путать со средним межкластерным расстоянием, которое мы тоже обозначали $F_1$ выше) будет просто средним гармоническим гомогенности и полноты:
$$ V_{1} = frac{2 cdot Homogeneity cdot Completeness}{Homogeneity + Completeness} $$
V-мера комбинирует в себе гомогенность и полноту таким образом, чтобы максимизация итоговой метрики не приводила к тривиальным решениям.
Коэффициент силуэта
Ещё одна метрика кластеризации, на этот раз уже не требующая разметки, это коэффициент силуэта (silhouette coefficient). Изначально коэффициент определяется для каждого объекта выборки, а метрика для результатов кластеризации всей выборки вводится как средний коэффициент силуэта для всех объектов выборки.
Чтобы ввести коэффициент силуэта $S(x_i)$, нам потребуются две вспомогательные величины. Первая, $A(x_i)$, — это среднее расстояние между $x_i$ и объектами того же кластера. Вторая, $B(x_i)$, — это среднее расстояние между $x_i$ и объектами следующего ближайшего кластера. Коэффициент силуэта вводится следующим образом:
$$ S(x_i) = frac{B(x_i) – A(x_i)}{max(B(x_i), A(x_i))} $$
В идеальном случае объекты «родного» кластера $x_i$ должны быть ближе к $x_i$, чем объекты соседнего кластера, то есть $A(x_i) < B(x_i)$. Однако это неравенство выполняется далеко не всегда. Если «родной» кластер $x_i$, например, имеет форму очень протяжённой ленты или просто большой размер, а недалеко от $x_i$ есть кластер поменьше, может оказаться, что $A(x_i) > B(x_i)$. Таким образом, если мы посмотрим на разность $B(x_i) – A(x_i)$, она может оказаться как положительной, так и отрицательной, но в идеальном сценарии всё же следует ожидать положительное значение. Сам же коэффициент $S(x_i)$ принимает значения от –1 до +1 и максимизируется, когда кластеры кучные и хорошо отделены друг от друга.
Коэффициент силуэта особенно полезен (по сравнению с другими приведёнными метриками) тем, что одновременно и не требует разметки, и позволяет подбирать количество кластеров. См. подробнее в примере из документации scikit-learn.
Различия и выбор метрик качества кластеризации
Подводя итог в теме метрик качества в задаче кластеризации, отметим, что есть несколько разных сценариев использования этих метрик. Если вы уже определились с количеством кластеров, можно оптимизировать среднее внутрикластерное или среднее межкластерное расстояние. Если у вас ещё и есть разметка — гомогенность и полноту. V-мера за счёт сочетания гомогенности и полноты в целом позволяет выполнять и подбор количества кластеров.
Однако разметка, с одной стороны, есть далеко не всегда, а с другой стороны, в задаче кластеризации часто очень субъективна. Сложность кластеризации в том, что на одной и той же выборке нас вполне могут устроить сразу несколько различных вариантов кластеризации, то есть задача поставлена некорректно и имеет более одного решения. Формализовать, какие решения нас устроят, на практике довольно сложно, поэтому сама по себе задача кластеризации решается не слишком хорошо.
Если разметки нет и число кластеров не фиксировано, лучшая метрика на практике — коэффициент силуэта. Исключение — ситуация, когда результат кластеризации используется далее для решения некоторой задачи обучения с учителем (как было в примере классификации изображений с помощью visual bag of words). В этом случае можно абстрагироваться от качества кластеризации и выбирать такой алгоритм кластеризации и такие его гиперпараметры, которые позволят лучше всего решить итоговую задачу.
Время на прочтение
5 мин
Количество просмотров 29K
Кластеризация — разбиение множества объектов на подмножества, называемые кластерами. Кластеризация, будучи математическим алгоритм имеет широкое применение во многих сферах: начиная с таких естественно научных областей как биология и физиология, и заканчивая маркетингом в социальных сетях и поисковой оптимизацией.
Существует множество алгоритмов кластеризации, однако ниже будет рассмотрен метод k-средних, так как он является наиболее лаконичным и простым для понимания.
Кластеризация методом k-средних:
Исходной задачей будет распределение произвольного количества n-мерных точек по k кластерам.
-
Случайным образом создаются k точек, в дальнейшем будем называть их центрами кластеров;
-
Для каждой точки ставится в соответствии ближайший к ней центр кластера;
-
Вычисляются средние арифметические точек, принадлежащих к определённому кластеру. Именно эти значения становятся новыми центрами кластеров;
-
Шаги 2 и 3 повторяются до тех пор, пока пересчёт центров кластеров будет приносить плоды. Как только высчитанные центры кластеров совпадут с предыдущими, алгоритм будет окончен.
Приступим к реализации алгоритма:
Исходные данные алгоритма:
-
n — количество строк;
-
k — количество кластеров;
-
dim — размерность точек (пространства).
Выходные данные алгоритма:
-
cluster — двумерный массив размерностью dim * k, содержащий k точек — центры кластеров;
-
cluster_content — массив, содержащий в себе k массивов — массивов точек принадлежащих соответствующему кластеру.
def clusterization(array, k):
n = len(array)
dim = len(array[0])
cluster = [[0 for i in range(dim)] for q in range(k)]
cluster_content = [[] for i in range(k)]
for i in range(dim):
for q in range(k):
cluster[q][i] = random.randint(0, max_cluster_value)
cluster_content = data_distribution(array, cluster)
Переменные заданы. Первичные центры кластеров созданы с помощью библиотеки random(стр. 9 – 11) max_claster_value — константа задающая примерные границы исходного множества;
При помощи функции data_ditribution() произведено первичное распределения точек по кластерам (стр. 13). Рассмотрим эту функцию подробнее:
def data_distribution(array, cluster):
cluster_content = [[] for i in range(k)]
for i in range(n):
min_distance = float('inf')
situable_cluster = -1
for j in range(k):
distance = 0
for q in range(dim):
distance += (array[i][q]-cluster[j][q])**2
distance = distance**(1/2)
if distance < min_distance:
min_distance = distance
situable_cluster = j
cluster_content[situable_cluster].append(array[i])
return cluster_content
Для каждой строчки (стр. 5) высчитывается расстояние до каждого центра кластеров. Здесь применяется стандартный алгоритм:
-
За начальное кратчайшее расстояние (min_distance) берётся несоизмеримо большое со значениями точек число;
-
Затем происходит вычисление расстояния до центра каждого кластера;
Вычисление расстояния между точкой и центром в n-мерном пространстве. -
Если вычисленное расстояние меньше минимального, то минимальное расстояние приравнивается к вычисленному и точка привязывается к этому кластеру (situable_cluster);
-
После обработки точки, в массив cluster_content в выбранный кластер (situable_cluster) кластер вкладывается значение точки.
Функция возвращает массив cluster_content.
В дальнейшем, как и полагается, данная последовательность действий обращается в цикл:
privious_cluster = copy.deepcopy(cluster)
while 1:
cluster = cluster_update(cluster, cluster_content, dim)
cluster_content = data_distribution(array, cluster)
if cluster == privious_cluster:
break
privious_cluster = copy.deepcopy(cluster)
Данный цикл целостным образом описывает шаг 4 из описании алгоритм k-средних (см. выше).
После распределения точек по центрам кластеров происходит перераспределение уже центров кластеров по привязанным к ним точкам (стр. 2). Рассмотрим функцию cluster_content() подробнее:
def cluster_update(cluster, cluster_content, dim):
k = len(cluster)
for i in range(k): #по i кластерам
for q in range(dim): #по q параметрам
updated_parameter = 0
for j in range(len(cluster_content[i])):
updated_parameter += cluster_content[i][j][q]
if len(cluster_content[i]) != 0:
updated_parameter = updated_parameter / len(cluster_content[i])
cluster[i][q] = updated_parameter
return cluster
Для каждого кластера, для каждого из n измерений вычисляется новое значения с помощью незамысловатого среднего арифметического: стр. 8-9 — складываются все значения; стр. 11-12 — сумма делится на количество точек в кластере; стр. 13 — кластер принимает обновлённое значение.
На данном месте алгоритм заканчивает свою работу. Полный алгоритм выглядит следующим образом:
def clusterization(array, k):
n = len(array)
dim = len(array[0])
cluster = [[0 for i in range(dim)] for q in range(k)]
cluster_content = [[] for i in range(k)]
for i in range(dim):
for q in range(k):
cluster[q][i] = random.randint(0, max_cluster_value)
cluster_content = data_distribution(array, cluster)
privious_cluster = copy.deepcopy(cluster)
while 1:
cluster = cluster_update(cluster, cluster_content, dim)
cluster_content = data_distribution(array, cluster)
if cluster == privious_cluster:
break
privious_cluster = copy.deepcopy(cluster)
Перейдём к визуализации:
Визуализируем результат алгоритма для 3-х и 2-х мерного исходных пространств. Воспользуемся библиотекой mathplotlib:
import matplotlib.pyplot as plt
from mpl_toolkits import mplot3d
import numpy as np
Визуализация для 2-х мерного пространства происходит следующим образом:
def visualisation_2d(cluster_content):
k = len(cluster_content)
plt.grid()
plt.xlabel("x")
plt.ylabel("y")
for i in range(k):
x_coordinates = []
y_coordinates = []
for q in range(len(cluster_content[i])):
x_coordinates.append(cluster_content[i][q][0])
y_coordinates.append(cluster_content[i][q][1])
plt.scatter(x_coordinates, y_coordinates)
plt.show()
Grid() — создание сетки. xlabel(), ylabel() — названия осей. Затем в массивы, соответствующие осям вкладываются значения точек. После такой операции для каждого кластера вызывается функция scatter() — разброс точек по плоскости. В конце вызывается функция отображения — show().
Аналогичным образом визуализируется результат алгоритма для 3-х мерного пространства:
def visualisation_3d(cluster_content):
ax = plt.axes(projection="3d")
plt.xlabel("x")
plt.ylabel("y")
k = len(cluster_content)
for i in range(k):
x_coordinates = []
y_coordinates = []
z_coordinates = []
for q in range(len(cluster_content[i])):
x_coordinates.append(cluster_content[i][q][0])
y_coordinates.append(cluster_content[i][q][1])
z_coordinates.append(cluster_content[i][q][2])
ax.scatter(x_coordinates, y_coordinates, z_coordinates)
plt.show()
Визуализация пошагово, с отображением центров: https://drive.google.com/drive/folders/10WLlqfw__Bqr3ERjzZyMjk3ZNOIPPad2?usp=sharing
Таким образом, мы выполнили необходимые и достаточные условия для анализа и реализации кластеризации методом k-средних.
Кластеризация данных – это метод обнаружения и визуализации групп связанных между собой предметов. Данный инструмент часто используется в приложениях, обрабатывающих большие объемы данных. Кластеризация – пример обучения без учителя. В отличие от нейронных сетей или деревьев решений, алгоритмам обучения без учителя не сообщаются правильные ответы. Их задача – обнаружить структуру в наборе данных, когда ни один элемент данных не является ответом.
Постановка задачи
Совсем недавно, компания Adobe, выпустила на iPhone програмку Adobe Color CC. До меня дошла эта новость и я понял, что пора объяснить, как же это делается!) Метод известный и описан во многих местах, но наше – всегда лучше))) Задача у нас простая, методом K-средних сделать K кластеров цветов картинки и выявив центры кластеров, посмотреть какому цвету он соответствует. Таким образом мы узнаем основные тона картинки 😉 Но сначала, немножко теории.
Иерархическая кластеризация
Методоово кластеризации бывает несколько, в основном, мы будем рассматривать метод k-средних, но есть так же и метод иерархической кластеризации. Алгоритм иерархической кластеризации строит иерархию групп, объединяя на каждом шаге две самые похожие группы. В начале каждая группа состоит из одного элемента, в данном случае – одного блога. На каждой итерации вычисляются попарные расстояния между группами, и группы, оказавшиеся самыми близкими, объединяются в новую группу. Так повторяется до тех пор, пока не останется всего одна группа.
Иерархическая кластеризация в действии |
Здесь схожесть элементов представлена их относительным расположением – чем элементы ближе друг к другу, тем более они схожи. В начале каждый элемент – это отдельная группа. На втором шаге два самых лизких элемента A и B объединены в новую группу, расположенную посередине между исходными. На третьем шаге эта новая группа объединяется с элементом C. Поскольку теперь два ближайших элемента – это D и E, то из них образуется новая группа. И на последнем шаге две оставшиеся группы объединяются в одну.
Результаты иерархической кластеризации обычно представляются в виде графа, который называется дендрограммой. На нем изображено, как из узлов формировалась иерархия.
Дендрограмма как способ визуализации |
На дендрограмме представлены не только ребра графа, показывающие, из каких элементов составлен каждый кластер, но и расстояния, говорящие о том, как далеко эти элементы отстояли друг от друга. Кластер AB гораздо ближе к составляющим его элементам A и B, чем кластер DE к элементам D и E. Изображение графа таким образом помогает понять, насколько схожи элементы, вошедшие в кластер. Эта характеристика называется теснотой (tightness) кластера.
Иерархическая кластеризация дает на выходе симпатичное дерево, но у этого метода есть два недостатка. Древовидное представление само по себе не разбивает данные на группы, для этого нужна дополнительная работа. Кроме того, алгоритм требует очень большого объема вычисле- ний. Поскольку необходимо вычислять попарные соотношения, а затем вычислять их заново после объединения элементов, то на больших наборах данных алгоритм работает медленно. И тут на помощь приходит другой метод!)
Кластеризация методом k-средних
Он в корне отличается от иерархической кластеризации, поскольку ему нужно заранее указать, сколько кластеров генерировать. Алгоритм определяет размеры кластеров исходя из структуры данных.
Кластеризация методом K-средних начинается с выбора k случайно расположенных центроидов (точек, представляющих центр кластера). Каждому элементу назначается ближайший центроид. После того как назначение выполнено, каждый центроид перемещается в точку, рассчитываемую как среднее по всем приписанным к нему элементам. Затем назначение выполняется снова. Эта процедура повторяется до тех пор, пока назначения не прекратят изменяться.
На первом шаге случайно размещены два центроида (изображены темными кружками). На втором шаге каждый элемент приписан к ближайшему центроиду. В данном случае A и B приписаны к верхнему центроиду, а C, D и E – к нижнему. На третьем шаге каждый центроид перемещен в точку, рассчитанную как среднее по всем приписанным к нему элементам. После пересчета назначений оказалось, что C теперь ближе к верхнему центроиду, а D и E – по-прежнему ближе к нижнему. Окончательный результат достигается, когда A, B и C помещены в один кластер, а D и E – в другой.
В данном методе, ОЧЕНЬ важен выбор центроидов, т.к. в конечном итоге все точки к ним стянутся, но если все 4 будут очень близкие цвета и из них будут состоять большие области, то в конечном итоге вы и получите приблизительно одинаковые центры. Таким образом центроиды должны быть максимально разнесены (желательно равномерно и равноудаленно).
Для большего понимания приведу пример. Давайте рассмотрим такую картинку.
Чтобы решить поставленную задачу в такой, монотонной четко разделенной картинке, достаточно написать небольшой код:
def k_cluster_by_area(points, k=4): clusters_by_area = {} for point in points: clusters_by_area[point] = clusters_by_area.setdefault(point, 0) + 1 sorted_x = sorted(clusters_by_area.items(), key=operator.itemgetter(1)) return sorted_x[-k:]
Тут мы просто посчитали, сколько раз, какой цвет встречается, отсортировали и вывели k результатов. Этот метод можно применить к любому изображению, работает достаточно быстро, но есть НО! А именно, мы не получим ТОНАЛЬНОСТЬ картинки! Однако, как я писал выше, в нашей задаче важны центры, вот эту задачу и решает этот код!) Мы берем самые популярные цвета, за изначальные центры кластеров, а потом они постепенно мутируют) На самом деле, желательно брать не просто первые k, а посчитать все комбанации расстояний и выбрать k по макс расстояниям, чтобы центры не оказались рядом (например 2 оттенка серого). Но мне влом и вы сделаете это сами)
Надеюсь, стало яснее. Таким образом, как найти изначальные центры кластеров нам известно, осталось только свести все пиксели в эти точки. А теперь сам код:
def k_cluster(img, distance=euclidean, k=4): points = get_points(img) clusters = {} # центры кластеры cluster_centers = k_cluster_by_area(points, k) for cluster_center in cluster_centers: # avg_sum_coords - список сумм r, g, b координат всех пикселей кластера # count_points - кол-во пикселей в кластере # какие конкретно пиксели в кластере не запоминаем в целях экономии памяти и времени clusters[cluster_center] = {'avg_sum_coords':list(cluster_center), 'count_points':1} while len(points): # находим ближайшие точки к центрам кластеров bestmatches = {} for point in points: # расстояние каждой точки до каждого центра кластера for cluster_center in clusters: d = distance(cluster_center, point) if not bestmatches: bestmatches = {'distance':d, 'point': point, 'cluster_center': cluster_center} if d < bestmatches['distance']: bestmatches = {'distance': d, 'point': point, 'cluster_center': cluster_center} # перекидываем ближайшие точки в соответствующий кластер bestmatch = bestmatches['point'] cluster_center = bestmatches['cluster_center'] points2 = [] new_points = [] for point in points: if point == bestmatch: new_points.append(point) else: points2.append(point) points = points2[:] # меняем центр кластера print clusters.keys(), len(points) r_avg, g_avg, b_avg = 0, 0, 0 cluster_avg_coord = clusters[cluster_center]['avg_sum_coords'] new_points_len = len(new_points) clusters[cluster_center]['count_points'] += new_points_len claster_count_points = clusters[cluster_center]['count_points'] cluster_avg_coord[0] += bestmatch[0]*new_points_len cluster_avg_coord[1] += bestmatch[1]*new_points_len cluster_avg_coord[2] += bestmatch[2]*new_points_len new_center = cluster_avg_coord[0]/claster_count_points, cluster_avg_coord[1]/claster_count_points, cluster_avg_coord[2]/claster_count_points # если новый центр кластера отличается от предыдущего if not new_center in clusters: clusters[new_center] = clusters.pop(cluster_center) return clusters.keys()
Работает это дело не быстро ;(, поэтому я советую сжимать картинку, хотя бы до разрешения 100 на 100.
Понятно, что идея не нова и сравнить вы можете, например с данным сайтом http://design-seeds.com
Давайте разберем пример. С указанного выше сайта взяли картинку и их разбиение на цвета. Здесь можно заметить, что их задача – это выбрать основные центровые цвета, что немного отличается от нашей задачи, но все же.
Пример с сайта |
Сначала посмотрим на наш метод выбора первоначальных центров, исходя из предположения простого количественного показателя использования цветов. Изначальные центры ниже подписано, сколько раз цвет встречался на картинке.
7 раз |
4 раза |
4 раза |
4 раза |
5 раз |
6 раз |
Ммммм, кто бы сказал по картинке, что это самые часто встречающиеся цвета? Понятно, что это вызвано так же и сжатием до 100 пикселей по меньшей из сторон, но все же. Сразу можно сказать, что притягивать к таким центрам остальные цвета сложно. В результате:
Все плохо, как и следовало ожидать. А что если попробовать представить спектр цветов в виде прямой от 0 до 255, равномерно разбить его на k отрезков и сделать центры точками с одинаковыми равномерными координатами? По понятным причинам они все будут оттенками серого. Пример для 6:
В результате:
Вот, уже лучше) Однако, чем плох этот метод? А тем, что этих цветов может и не быть на картинке (она может просто красная), получается, что изначально мы как бы сбиваем мишень(
Как вывод, могу сказать, что метод k-средних очень чувствителен к этим самы центроидам. И их поиск – это не менее важная задача, которую я Вам и предлагаю решить) Например, представьте что вы как бы объединяете точки, расширяя территории. Объединяя их в 3х мерном пространстве линейно увеличивая координаты. И без учета количественных показателей весов!)