Алгоритм минимального покрывающего дерева (MST) [Кластерный анализ и Python]

В продолжение цикла заметок о кластерном анализе сегодня поговорим о алгоритме минимального покрывающего дерева (MST) и работе с ним на языке python.



В первой заметке уже был описан набор данных, который мы исследуем, а также приведён список литературы на тему, поэтому перейдём сразу к делу.

Теория

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

Алгоритм минимального покрывающего дерева (MST) строит граф из N-1 рёбер так, чтобы они соединяли все N точек и обладали минимальной суммарной длиной. Такой граф называется кратчайшим незамкнутым путём, минимальным покрывающим деревом или каркасом графа.

Описание работы алгоритма:

• Найти пару точек с наименьшим расстоянием (весом ребра) и соединить их ребром;

• пока в выборке остаются изолированные точки:

• найти изолированную точку, ближайшую к некоторой неизолированной;

• соединить эти две точки ребром;

• удалить K-1 самых длинных рѐбер;

На последнем шаге алгоритм удаляет из полученного графа K-1 самых длинных рёбер, что приводит к распаду графа на отдельные связанные компоненты, которые и будут искомыми кластерами.

Общеизвестны два недостатка этого алгоритма:

• ограниченная применимость. Алгоритм наиболее подходит для выделения кластеров типа сгущений или лент. Наличие разреженного фона или «узких перемычек» между кластерами приводит к неадекватным результатам;

• высокая трудоёмкость — для построения кратчайшего незамкнутого пути требуется O(N^3) операций.

Практика

Для работы с графами будет использоваться NetworkX — библиотека для создания и манипуляции графами. Кроме этого, она предоставляет готовую реализацию алгоритма нахождения минимального остовного дерева и удобное API для визуализации результатов.

Для начала необходимо получить исходные данные и расстояния между всеми элементами набора.

from scipy.spatial.distance import pdist
names, data = get_data()
dist = pdist(data, 'euclidean')

Функции get_data() и pdist() были описаны в первой заметке, поэтому не будем на них останавливаться.

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

    from collections import deque
    import networkx as nx 
    s = nx.Graph()
    s.add_nodes_from(names)
    dq = deque(dist)
    len_x = len(names)
    for x in xrange(len_x - 1):
        for y in xrange(x + 1, len_x):
            s.add_edge(names[x], names[y], weight=dq.popleft())

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

    import matplotlib.pyplot as plt
    mst = nx.minimum_spanning_tree(s)
    plt.hist([edge[2]['weight'] for edge in mst.edges_iter(data=True)], 100, color='red', alpha=0.3) 

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

    r = nx.Graph()
    r.add_nodes_from(names)
    edges = [edge for edge in mst.edges_iter(data=True) if edge[2]['weight'] <= 0.05]
    r.add_edges_from(edges)

Как вы могли заметить, в отличии от классического алгоритма MST, я предпочёл отрезать рёбра не по количеству кластеров, а по пороговому уровню, как в иерархическом алгоритме. Это позволило снизить влияние на результат кластеризации качество предварительной подготовки данных.

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

    def graph_draw(g):
        """Draw graph"""
        plt.figure()
        nx.draw_graphviz(g, with_labels=False, node_size=3, prog='neato') 
    graph_draw(mst)
    graph_draw(r)
    plt.show()

Минимальное остовное дерево графа

Результат кластеризации