НОВОСТИ Агломеративная кластеризация: алгоритм, быстродействие, код на GitHub

Alvaros
Онлайн
Регистрация
14.05.16
Сообщения
21.452
Реакции
101
Репутация
204



Несколько лет назад мне потребовалось очень качественно кластеризовать относительно неплотные графы среднего размера (сотни тысяч объектов, сотни миллионов связей). Тогда оказалось, что алгоритма с подходящим набором свойств просто не существует, несмотря на всё разнообразие методов, придуманных человечеством за многие десятилетия. Имеющиеся решения работали либо просто очень плохо, либо очень плохо и к тому же медленно.


К счастью, оказалось, что идеи, заложенные в метрики качества кластеризации, о которых я , можно адаптировать и создать на их основе алгоритм кластеризации. Он достигает очень высоких показателей качества и к тому же работает очень быстро благодаря некоторым удачным аналитическим свойствам оптимизируемых величин. Алгоритм относится к классу агломеративных: основной операцией является слияние нескольких уже имеющихся кластеров в один более крупный кластер.


Об этом алгоритме и пойдёт речь в статье. Под катом читателей ждут математическое описание алгоритма, техники уменьшения его временной сложности, код на GitHub и модельные наборы данных.


1. Агломеративная кластеризация



Пусть даны множество объектов
814e648700be7248655d674954ec4b32.svg
и функция попарной близости объектов
d7705a68b5e415c9438b71d5425656d5.svg
. Каждый элемент близок сам себе, а функция близости симметрична:



8a8db57fc1f44f3a7b79e593cdf3aebb.svg




c8a6d01635a27022667ade7fc26ffbfe.svg



Также определено некоторое разбиение множества
493c1c008018df9bed4910321f29ff00.svg
на кластеры:



c4b199f11417dcfccbd4966badeb84b4.svg



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


Множество кластеров будем обозначать
028e8a1e74b8335d284dc41e5a41cf70.svg
:



64552944475b76bd9fa6247b1363e846.svg



Пусть выбраны объект
c6a200889233190d0ffe86be0f9fadc8.svg
и множество
d11f6cf94553babb595dfc5abd9b06ba.svg
. Воспользуемся интуицией BCubed-метрик качества кластеризации и определим для этой пары значения точности и полноты:



088e28a4d957a1218842be9236a6ebc1.svg




b1e387719fd88174fa41e3225484c1f1.svg



Заметим, что эти определения удобны тем, что
817b92407f764f57af9226e50cc788fd.svg
не обязательно находится во множестве
6d6a4f78fbacd6edecc018ce8ad3e364.svg
.


Теперь можно вычислить и средние значения точности и полноты для элементов
6d6a4f78fbacd6edecc018ce8ad3e364.svg
относительно самого
6d6a4f78fbacd6edecc018ce8ad3e364.svg
:



ab203f755e220e2527afdfe477671f4c.svg




96b41a4260370be49fcf697cb45f5659.svg



Определим композицию этих показателей, воспользовавшись идеей метрики ECC:



14c337915e2b3d0104e2bbd6d63d8ffe.svg



Здесь параметр
7234a52ba041cdb09b9328a047048fb2.svg
, как и в
3c9bcade02ed1bcb586c485ae9c4eba8.svg
-меры, задаёт относительную важность полноты относительно точности. В простейшем случае
15903772463609f1294d9587d75455bb.svg
и показатели равноправны.


Рассмотрим два непересекающихся подмножества
e97eec3dab7c27b0b911e9ce196f4e7a.svg
и
c30fd554ca506c6c2a613107a5143dc0.svg
. Если они являются кластерами, то средние точность и полнота входящих в них элементов равны:



bd98b3bc9a250f94cb976bf9eb60f6fe.svg




0e45709c6e2b614f922fc93144106ed7.svg



Определим и композицию метрик для этого случая:



9777f12be81c944852b8b344270d3dfa.svg



Если объединить элементы множеств
6d6a4f78fbacd6edecc018ce8ad3e364.svg
и
c62ff25ef4caeaeaef7122a489ef9d07.svg
в один кластер, показатели станут равны:



094ae883662c5dcf450a7e2304b4b5a0.svg




26d0e7040d12aca24084e27ce20487ab.svg



А композиция метрик определится просто как
be0c3b4bc5c657cb11fb85971bf82c37.svg
.


При объединении кластеров
6d6a4f78fbacd6edecc018ce8ad3e364.svg
и
c62ff25ef4caeaeaef7122a489ef9d07.svg
композиция метрик для входящих в них элементов изменится на величину



8ef479bc0f2e704ce2c8393876cc2498.svg



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

Даны:
  • Множество объектов
    ca0945c57d7f91df83f74b40edb1c618.svg
  • Симметричная попарной близости
    d7705a68b5e415c9438b71d5425656d5.svg
  • Тривиальное множество кластеров:
    72715edda7d4822ac501be47b231f0e6.svg
    ,
    375bd40c94d9d4929350e5bac1e8f3e3.svg


Основной цикл:
  1. Если
    4d46d6689a3bbd8c8cd838d777322e76.svg
    , вернуть
    028e8a1e74b8335d284dc41e5a41cf70.svg
    и завершить алгоритм.
  2. Выбрать
    a6f4e8b5efc92346cbe1a8fa761b438b.svg
    .
  3. Если
    985741561f27d13cee9a6facfb4fbac2.svg
    , вернуть
    028e8a1e74b8335d284dc41e5a41cf70.svg
    и завершить алгоритм.
  4. В противном случае положить
    dc2948e520a0e79c1fa205e6404ee330.svg
    и вернуться к шагу 1.

Основной проблемой такого алгоритма будет быстродействие. Каждая итерация основного цикла требует вычисления функции
c14a59bde821353b1d40d63f5578b34d.svg
для всех пар кластеров. Но оказывается, что эту часть алгоритма можно существенно ускорить.

2. Ускорение алгоритма



Сложность алгоритма проявляется в двух аспектах. Первый связан с тем, что количество пар кластеров на шаге 2 является квадратичным от общего числа имеющихся кластеров; второй — с тем, что вычисление величины
d1c1387edd13f82aba0e19f2479c91c2.svg
при наивной реализации требует времени, пропорционального как минимум
234c85b1c091003a644d48c588ccf173.svg
.


Ускорить алгоритм в связи с первым аспектом достаточно просто. Заметим, что после объединения кластеров
be9fcb960a47cc7659becd7ce4c851d5.svg
и
c1355d4ee682504c5920c8c5cc9bd51b.svg
значения функции
c14a59bde821353b1d40d63f5578b34d.svg
нужно обновить лишь на тех парах кластеров, в которые входит либо
be9fcb960a47cc7659becd7ce4c851d5.svg
, либо
c1355d4ee682504c5920c8c5cc9bd51b.svg
.


Другими словами, величина
6eb600875853bee46568f0198c2c8f0b.svg
не меняется после объединения кластеров
be9fcb960a47cc7659becd7ce4c851d5.svg
и
c1355d4ee682504c5920c8c5cc9bd51b.svg
, если
cfcb76f862e62efbb6c2aaec388eed0f.svg
. Следовательно, обновится лишь линейное, а не квадратичное по размеру
028e8a1e74b8335d284dc41e5a41cf70.svg
число значений.


Хранить их при этом можно в любой set-подобной структуре, так что извлечение максимального элемента будет занимать
edac6897ca5bab63e49737e4cc6421e3.svg
времени.


Справиться со вторым аспектом несколько сложнее. Чтобы лучше понять, как здесь работает ускорение, рассмотрим произвольную аддитивную по ячейкам матрицы функцию
12ab92a586fb3d91b2b9c1e1ca0dc58c.svg
. Пример такой функции — простая сумма значений.


Пусть множество
493c1c008018df9bed4910321f29ff00.svg
разбито на три кластера. Тогда функция
12ab92a586fb3d91b2b9c1e1ca0dc58c.svg
определена на всех девяти прямоугольных фрагментах матрицы, которые соответствуют парам различных кластеров. Нас интересует, как обновить её значения при объединении двух кластеров в один.


zig3buip1wktyefhbswiqnf4iak.png



После такой операции останется четыре пары кластеров и четыре интересующих нас значения функции. В силу её аддитивности эти четыре значения вычисляются из предшествующих за константное время, не зависящее от размеров кластеров. Правила такого пересчёта показаны на картинке выше.


Теперь заметим: если в ячейках матрицы находятся близости элементов, а функция
12ab92a586fb3d91b2b9c1e1ca0dc58c.svg
осуществляет простое суммирование, мы получаем правила для вычисления и обновления характеристик точности. Средние показатели метрики точности для элементов из
40aa32341abe8afc6e84b17d8dc6f99b.svg
до и после объединения будут равны, соответственно:



c2fb86f3ed96410d2018bc4dfab6fded.svg




24ebce325cb28d53c9027322d4c6900e.svg



Аналогично можно поступить и с полнотой. Определим аддитивную функцию
bdba7499baaa2899811e34409321d6eb.svg
так, что она будет суммировать нормированные по строкам близости из исходной матрицы. Таким образом:



72527c348f6b86bca2870c1cac2d3109.svg



Сумма значений
bdba7499baaa2899811e34409321d6eb.svg
по элементам одной строки матрицы всегда равняется единице. Тогда:



eb8214e7c281863557e603c257d19fe9.svg




90059c5042c384e3f22b3f6cc6ed4220.svg



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


Аналогично можно вычислить показатели, необходимые для определения величины
fd7922a785a77a49b96b001df5ea1920.svg
:



939b5f4bc4a199eb14918d046578ebbb.svg




3a11d558fab162533be54f0cb4173b7d.svg




8ddf4a25eccc5d915d230f469bec193a.svg




8af3e3cda4ac33605096ca27d8bc3573.svg



То есть благодаря аддитивности каждой введённой величины по элементам матрицы все вычисления работают за константное время! Обновление всех значений функции в худшем случае потребует
a8c3d4a6e963aa33224286ee2cfdfe56.svg
времени.


Поскольку каждая итерация уменьшает число кластеров на единицу, полное выполнение алгоритма потребует
5a6905e58424db44945e46f3cb1f367a.svg
времени. Если в среднем по всем итерациям один кластер оказывается связан с
16da507b2fc389688ef0659939dcc647.svg
другими, оценка временной сложности улучшается до
c19b0d2ff5b112877f306646f1d70294.svg
.


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

4. Реализация



Описанный алгоритм реализован на C++ и выложен под лицензией MIT . Реализация компилируется и запускается на Linux и Windows.


Программу легко собрать:


git clone
cmake .
cmake --build .


После этого её можно запустить на одном из примеров, находящихся в том же проекте:


./agglomerative_clustering < data/2d_similarities.tsv > 2d_clusters.tsv


В результате будет сгенерирован файл, формат которого совпадает с форматом, в котором описывается разметка кластеризации (см. )

5. Наборы данных



Первый поставляемый с программой набор данных — синтетический. 18 343 точки на плоскости получены из 1000 нормальных двумерных распределений со случайными целочисленными центрами. Близость между двумя точками определяется по формуле:



cd2a6ad4dc550b77aef38cc3b85e8af0.svg



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


ze6khfakhukhp65qpfaht2plvjm.png

[SUP][SUB]Синтетический набор точек на плоскости[/SUB][/SUP]


Алгоритм кластеризации должен быть устойчив к потере информации о связях между объектами. Чтобы продемонстрировать это свойство, в наборе из всех близостей, существенно превосходящих ноль, была оставлена примерно одна треть (228 018 штук). Увеличивая параметр -f (который соответствует параметру
7234a52ba041cdb09b9328a047048fb2.svg
в этой статье), можно добиться практически тех же показателей качества, что и на графе до удаления рёбер:


./agglomerative_clustering -f 10 < data/2d_similarities.tsv > 2d_clusters.tsv
../cluster_metrics/cluster_metrics data/2d_markup.tsv 2d_clusters.tsv
ECC 0.62411 (0.62411)
BCP 0.74112 (0.74112)
BCR 0.68095 (0.68095)
BCF1 0.70976 (0.70976)


Здесь для вычисления метрик используется программа .


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


Чтобы отчасти объяснить, как могут выглядеть данные, с которыми приходится встречаться в реальных приложениях, я собрал . Он получен в ходе решения одной из аналитических задач, предоставляется без разметки и служит исключительно для проверки быстродействия алгоритмов.


cmake -DCMAKE_BUILD_TYPE=release .
cmake --build .
tar zxvf data/doc.similarities.tar.gz
./agglomerative_clustering < doc.similarities.tsv > doc.clusters.tsv


Программа в процессе работы выводит информацию о времени выполнения разных этапов:


loaded 5000000 pairs
loading documents: finished in 11.323s
clustering documents: finished in 28.653s
printing clusters: finished in 0.038s


В этом примере задача кластеризации 301 836 элементов с пятью миллионами связей занимает около полуминуты на моей виртуальной машине. Такого быстродействия более чем достаточно в приложениях, для которых разрабатывался алгоритм. Программу можно дополнительно ускорить за счёт параллелизма и хардкорных оптимизаций, но рассказ о них тянет на отдельный пост.
 
Сверху Снизу