"Видеоурок"

Свежий взгляд на алгоритмы кластеризации

Глубокое погружение в новую методологию распознавания кластеров.

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

Постановка проблемы

Позвольте мне начать с небольшого контекста. Поскольку Medium - это платформа для блогов, а не место для научных исследований, я могу позволить себе рассказать историю, а не представлять результаты в очень сухой форме. Пожалуйста, простите меня, если это вас разочаровывает - это, по сути, своего рода «шоу и рассказ».

Пару лет назад я пытался решить задачу анализа растровых файлов. Если вы никогда раньше не работали с ними, посмотрите мою серию «Геопространственные приключения» для более подробной информации:



Геопространственные приключения. Шаг 5. Покинуть равнину или полетать над морем многоугольников.
Как следует из названия, пора нам сделать что-то серьезное и начать летать с помощью «Лазеров. Буквально. Мы… в направленииdatascience.com »



В двух словах - это матрицы чисел высоты, полученные в результате полета над областью, которую вы пытаетесь анализировать, с регулярным направлением на нее лазерного луча и измерения отклика. Полученную матрицу можно использовать для извлечения информации о различных объектах на поверхности. Эти объекты могут быть чем угодно, от гор и деревьев до зданий и мостов. Я был особенно заинтересован в извлечении информации о зданиях, а это означало, что у меня были коллекции плотно упакованных точек, которые мне почему-то нужно было интерпретировать как многоугольники. Если у вас есть одинокое отдельно стоящее здание посреди поля - это просто, все, что вам нужно сделать, это применить алгоритм выпуклой оболочки из коробки (из библиотеки shapely. Опять же, если вы новичок в этом - взгляните у моего букваря здесь). Но что, если мы смотрим на густонаселенную центральную часть города?

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

Теперь, когда мы столкнулись с кризисом COVID и у меня было много свободного времени во время поиска работы, я подумал, что было бы неплохо вернуться к проблеме и обобщить алгоритм. Я представляю здесь результат.

Уф… Это было довольно долгое вступление. Давайте погрузимся.

Одна из фундаментальных проблем кластеризации заключается в том, что мы действительно хорошо оцениваем ее визуально и редко задумываемся о ее механике. В результате действительно легко найти способы сбить с толку случайного наблюдателя и заставить его сделать ложный вывод. Если я увеличиваю набор точек на 2D-графике - я могу сказать, что они справедливо распределены, если я уменьшаю масштаб - они появляются в виде кластера. Итак, с самого начала нам нужно определить масштаб проблемы, это будет одним из основных параметров модели.

Самый простой способ сделать это - определить размер элементарной ячейки в пространстве нашей задачи. Если это простое 2D-изображение, это означает, что мы распределяем все точки в нашем наборе данных в ограничивающем квадрате, содержащем этот набор данных, и рисуем квадратные ячейки вокруг каждой точки. Если мы говорим о более сложной n-мерной проблеме, например, о наборе текстовых представлений word2vec - это будет ячейка единичного объема в n-мерном пространстве. Гораздо сложнее выявить группы данных, когда вы не видите их графического представления, поэтому действительно важно иметь правильный способ определения шкалы.
Если у вас есть категориальные переменные, их всегда можно сократить до набора двоичные переменные (1/0), или вы можете выбрать более сложную версию с использованием действительных чисел, если хотите построить между ними какие-то значимые отношения. В любом случае приблизительное количество требуемых ячеек / делений для каждого измерения может быть построено как:

Где N - количество точек (векторов) в наборе данных, n - количество измерений (непрерывные переменные), а k - количество категориальных переменных. . Я предполагаю, что категориальные переменные уменьшены до 1/0.

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

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

Алгоритм

Чтобы дать вам очень сжатую, краткую версию того, что мы делаем:

  1. Нормализуйте данные, чтобы каждая переменная находилась в интервале от 0 до 1.
  2. Умножьте каждую переменную на D и возьмите целую часть. Это эффективно стандартизирует каждую точку, помещая ее в определенную ячейку сетки.
  3. Подсчитайте количество точек в каждой ячейке сетки и отбросьте ячейки с номерами ниже определенного порога. Это еще один важный параметр, который мы собираемся передать в модель. В этой реализации я смотрю на конкретные процентили, но им можно управлять с помощью произвольного числа.
  4. Результирующая структура может быть представлена ​​в виде набора деревьев (или, скорее, графиков, но давайте кататься с деревьями), и все, что нам нужно сделать для идентификации наших кластеров, - это запустить алгоритм поиска в ширину (BFS) для подсчета экземпляров несвязанных деревьев.
  5. Наконец, как только основные кластеры были идентифицированы, необходимо отнести точки, которые не были учтены, поскольку их ячейки не прошли пороговое значение. Здесь я делаю это вне основной функции, вычисляя минимальное расстояние до центроидов идентифицированных кластеров. Это не единственный вариант, и я считаю, что даже не лучший вариант. Но пока это подойдет, и (!) Его легко и быстро реализовать.

Хорошо, приступим к делу.

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

Загрузка необходимых нам библиотек:

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
pd.plotting.backend='seaborn'
import time
import itertools

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

Matplotlib и Seaborn для визуализации, также использующие Seaborn в качестве бэкэнда для pandas. Использование модуля времени для просмотра сравнительной производительности и использование itertools для предотвращения циклов (ну, по крайней мере, в явной форме Python).

Набор данных:

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

Мы собираемся их нормализовать:

А теперь вперед и визуализируйте:

plt.figure(figsize=(10, 10))
sns.scatterplot(ttf.x, ttf.y, alpha=0.5)

Милая.

В нашем случае количество ячеек можно рассчитать как:

D = int(np.sqrt(len(ttf_norm)))
tolerance = 50

Это дает нам D = 63. Таким образом, параметры модели, которые мы можем использовать, будут D = 60 (просто хорошее круглое число), и мы будем искать 50-й процентиль с точки зрения количества точек на ячейку (допуск). В целом - выбор толерантности зависит от того, насколько хорошо распределены наши данные. Если нет возможности проверить это визуально или иным образом - неплохо запустить модель с разными уровнями допуска и посмотреть, сколько кластеров найдено для разных уровней и какой результат имеет наибольший смысл (является наиболее стабильным).

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

X = ttdf_norm[['x', 'y']].values

Хорошо, мы готовы к строительству. Я собираюсь пройти через это шаг за шагом.

  1. Масштаб и дискретность:
Xint = (X * D).astype(int)

2. Создайте идентификаторы строк для каждой ячейки:

Xstr = np.apply_along_axis(
        ''.join,
        1,
        np.char.zfill(
            Xint.astype(str),
            len(str(D))
        )
    )

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

3. Создайте словарь с подсчетом частоты появления для всех ячеек. Так же легко, как :

unique, counts = np.unique(Xstr, return_counts=True)
tree = dict(zip(unique, counts))

4. Рассчитайте порог:

mincount = int(
        np.percentile(
            np.array(
                list(
                    tree.values()
                )
            )[
                np.array(
                    list(tree.values())
                ) >1
            ],
            tolerance
        )
    )

Этот немного занят, так что давайте посмотрим, что происходит. Мы создаем массив numpy из всех значений в нашем дереве и отбрасываем все с одной точкой или меньше. Мы не хотим, чтобы они слишком сильно влияли на процентили, поскольку они не приносят никакой ценности. Из полученного массива мы вычисляем значение, соответствующее процентилю, установленному нашей переменной допуска. В простейшем случае 50% соответствует медиане.

5. Обрезка дерева: избавьтесь от всего, что ниже mincount. Использование простого понимания dict:

tree = {k: v for k, v in tree.items() if v > mincount}

В качестве побочного примечания - при игре с numba.jit мне пришлось заменить его простым for loop на Python (мое сердце истекло кровью, когда я это сделал) - поскольку numba не может обрабатывать понимание списка / словаря в текущей реализации. Теперь, подумав об этом - это можно сделать вообще без циклов, используя массивы numpy для ключей и значений. Это должно ускорить код, хотя это не является основным узким местом, поэтому время расчета существенно не изменится.

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

По сути, нам нужно выполнить поиск в ширину. Это означает, что мы собираемся выбрать случайный идентификатор из нашего дерева, пометить его как посещенный и назначить всех его соседей в список 'tocheck' (очередь в классической реализации), затем мы выберем любое значение из списка 'tocheck' и повторим тот же процесс. В классической реализации это должно быть первое значение от «tocheck», следовательно, используются очереди, но в нашем случае нам придется пройти через все из них, поэтому порядок не имеет значения. Когда у нас заканчиваются элементы в списке «для проверки», мы завершаем наше первое дерево. Мы можем назначить его первому кластеру и выбрать следующий идентификатор из исходного списка.

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

Сначала мы генерируем список сдвигов для наших координат. Мы ищем всех ближайших соседей на расстоянии одного места. В случае 2D - это 9 точек (включая начало координат, в 3D - это 27 и т. Д. Мы должны иметь возможность выполнять сдвиги в любом количестве измерений, определяемом X. Здесь пригодится itertools.

l = list(itertools.product([-1, 0, 1], repeat=Xint.shape[1]))
l=list(map(list, l))
l.remove([0] * Xint.shape[1])

Здесь происходит следующее: сначала мы создаем все возможные перестановки -1, 0 и 1 в горизонтальном измерении матрицы X (в нашем случае 2, поскольку у нас есть только 2 переменные).

Они возвращаются в виде кортежей, поэтому нам нужно преобразовать их в списки - это задача второй строки.

Финальный шаг - убираем нулевые смещения, нам это не нужно.

Вот результат:

[[-1, -1], [-1, 0], [-1, 1], [0, -1], [0, 1], [1, -1], [1, 0], [1, 1]]

Теперь пусть lD будет длиной строкового представления нашего параметра D и ids быть единственным идентификатором, для которого мы пытаемся найти соседей. Я знаю, что использование id было бы более логичным, но у Python есть проблемы, поэтому вместо этого нам нужно использовать ids.

nbrs = np.array(
            [
                list(
                    [
                        int(
                            ''.join(
                                list(ids)[i * lD: (i + 1) * lD]
                            )
                        ) for i in range(len(l[0]))
                    ]
                )
            ] * len(l)
        ) + np.array(l)

Хорошо, здесь много чего происходит, так что разберемся: взятие списка строки - преобразование ее в список составляющих ее символов. Затем мы берем части этого списка с количеством элементов, соответствующим длине нашего параметра D (lD), и склеиваем их в строки с помощью соединения. Как только это будет сделано, мы снова преобразуем их в целые значения, и это дает нам копию соответствующей строки из матрицы X, то есть все переменные, соответствующие этой элементарной ячейке. Затем мы реплицируем этот список столько раз, сколько соседей мы пытаемся добраться (количество напряжений в списке l). Последний шаг - преобразовать итоговый список и список напряжений в массив numpy и сложить их. Это дает нам поэлементную сумму и, следовательно, координаты всех соседей. Обратите внимание, что это будет работать для любого количества переменных, поэтому это хорошо для n-мерных векторов. Проблемы потенциально начнутся, когда попадут в действительно большое количество переменных, но давайте пока не будем об этом беспокоиться.

Наконец, повторяем преобразование в идентификаторы:

nbrs = np.apply_along_axis(
            ''.join,
            1,
            np.char.zfill(
                nbrs.astype(str),
                lD
            )
        )

И вуаля, у нас есть идентификаторы соседей. Все, что осталось, - это быстрая проверка, что мы не посещали их раньше и действительно ли они принадлежат идентифицированному дереву.

Слушай, мама, никаких «петель» !!!

Эээ ... Вообще-то там было понимание одного списка, хотя оно довольно короткое. Тем не менее, если вы найдете способ обойтись без него - дайте мне знать!

Как только мы закончили подсчет деревьев, мы более или менее закончили. Все, что нам нужно сделать, это приписать единичные ячейки, которые мы не учли, а затем пометить соответствующие точки.

Давайте посмотрим на картинки!

Пришло время сравнить результаты выполнения нашей функции со стандартными методами кластеризации. Если вы хотите просмотреть последний, то вот очень хороший пост с целой кучей - 10!

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

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

Загрузите модель:

import sklearn
from sklearn.cluster import AffinityPropagation
model = AffinityPropagation(damping=0.9)

И запустите подгонку:

%%time
model.fit(X)

Первая строка - время выполнения.

Заключительный этап

ttf_norm['Affinity'] = model.predict(X)

На моем ноутбуке этот метод занял больше всего времени - около 21 секунды.

Тот же процесс с другими. Я немного поиграл с параметрами, чтобы получить наилучший результат прямо из коробки. В случае с Affinity - это не особо помогло.

Я должен подчеркнуть, что некоторые из методов - KMeans (прямой и мини-пакет), Spectral, Gaussian mix, Agglomerative и Birch - все имеют количество кластеров в качестве входных данных, что в некотором смысле противоречит цели. Вы должны знать количество кластеров перед их применением. На самом деле, мы можем сразу отказаться от них, поскольку они не пригодятся ни для одной интересной задачи, которую мы собираемся решить. Но давайте оставим их в покое, хотя бы для упражнения.

Заговор:

Вы можете видеть, что MeanShift, Birch, DBSCAN, Affinity работают плохо. Между прочим, DBSCAN оказался особенно сложным с точки зрения получения правильных параметров даже для перехода на этот этап. Смесь Гаусса выглядит лучшей из всех. Наш метод DenseTrees работает наравне с KMeans. Если мы изменим нашу атрибуцию второго пути, мы, вероятно, сможем получить такой же хороший результат, как и от GM (фактически, мы можем использовать GM для окончательной атрибуции). Между прочим, Gaussian Mixture была примерно в 6 раз быстрее, чем DenseTrees, так что и здесь есть возможности для улучшения (у нас все еще все в порядке с 300 мс - 6-е место). Главный результат здесь - мы превзошли все алгоритмы, у которых не было количества кластеров в качестве входных данных, и мы недалеко от лучших алгоритмов из тех, что имели. Не слишком потрепанный.

Заключительное примечание

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

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

Особые похвалы тем, кто может понять, где это.

Это просто точки, нанесенные в нормализованном масштабе. Все, что ниже 10 метров, обнуляется. Вы можете видеть, что мы поймали несколько деревьев в нашу сеть. К сожалению, мы мало что можем с этим поделать. По крайней мере на данный момент. Этот квадрат имеет разрешение 2000 на 2000 (всего 4 миллиона точек данных), и каждый пиксель соответствует квадрату 0,5 на 0,5 м. Чтобы уловить зазоры между зданиями, нам нужно работать с большой детализацией, поэтому для обработки я использовал D = 500 и допуск 10%. Обратите внимание, что критическое значение для этого было бы слишком большим, поскольку все точки уже расположены равномерно. На моем ноутбуке это заняло около 5 часов. Не хорошо. Исходный алгоритм, оптимизированный для 2D, занял около получаса, так что, безусловно, есть много возможностей для улучшения.

Алгоритм возвращает 247 кластеров. И вот они все во всей своей радужной красе:

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

На сегодня все. Жду вашей критики и предложений.