Участник:LKruglov/Алгоритм устойчивой кластеризации с использованием связей

Алгоритм устойчивой кластеризации с использованием связей
Последовательный алгоритм
Последовательная сложность	[math]\displaystyle{ O(n^{2,37} + n^2\log n) }[/math], если для вычисления количества связей используется алгоритм Копперсмита—Винограда
Объём входных данных	Вектор из [math]\displaystyle{ n }[/math] точек в виде векторов характеристик произвольного вида, дополнительно одно вещественное и одно челое число [math]\displaystyle{ k }[/math]
Объём выходных данных	Сгруппированные в [math]\displaystyle{ k }[/math] кластеров [math]\displaystyle{ n }[/math] классифицируемых точек представленных индексами во входном векторе

Авторы: Круглов Л. В.(1.1,1.2,1.3,1.6,1.8,1.10), Баранов М. С.(1.1,1.2,1.3,1.4,1.5,1.7,1.9,2.7)

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

Свойства и структура алгоритма

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

Алгоритм устойчивой кластеризации с использованием связей (robust clustering using links, ROCK [1]) построен на основе понятий связей между кластеризуемыми точками и их соседства. Две точки считаются соседями, если являются достаточно близким (схожими): [math]\displaystyle{ sim({p_i}, {p_j}) \gt = \theta }[/math], где [math]\displaystyle{ \theta }[/math] - заданное пороговое значение в интервале [math]\displaystyle{ (0, 1) }[/math]. Связь [math]\displaystyle{ link({p_i}, {p_j}) }[/math] между двумя точками определяется как число общих соседей между точками [math]\displaystyle{ {p_i} }[/math] и [math]\displaystyle{ {p_j} }[/math]. Таким образом, чем большее значение функции связи для двух точек, тем более вероятно они принадлежат одному кластеру.

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

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

Математическое описание алгоритма

Задача кластеризации сводится к выделению такого разбиения точек на кластеры, которое максимизирует сумму связей точек принадлежащих одному кластеру, и минимизирует эту функцию для точек из разных кластеров. Таким образом, для разбияния на [math]\displaystyle{ k }[/math] кластеров необходима максимизация следующей функции.

[math]\displaystyle{ E_l = \sum_{i=1}^k n_i * \sum_{p_q, p_r \in C_i} \frac{link(p_q, p_r)}{n_i^{1+2f(\theta)}} \qquad(1) }[/math]

Данная функция отличается от простой суммы связей для каждого кластера, так как не обеспечивает распределение точек, имеющих мало связей между собой, по разным кластерам. Поэтому действительную сумму связей в кластере следует разделить на ожидаемую сумму связей. Функция [math]\displaystyle{ f(\theta) }[/math] задается пользователем исходя из свойств обрабатываемых данных и ожидаемого характера кластеризации. Она должна обладать следующим свойством: каждая точка, принадлежащая кластеру [math]\displaystyle{ C_i }[/math], имеет приблизительно [math]\displaystyle{ n_i^{f(\theta)} }[/math] соседей в этом кластере.

Связь между двумя кластерами определяется как

[math]\displaystyle{ link({C_i}, {C_j}) = \sum^{}_{{p_q}\in{C_i}, {p_r}\in{C_j}}{link[{C_i}, {C_j}]} }[/math]

Тогда для выбора кластеров для объединения определяется целевая функция качества:

[math]\displaystyle{ g(C_i, C_j) = \frac{link[C_i, C_j]}{(n_i+n_j)^{1+2f(\theta)} - f_i^{1+2f(\theta)} - f_j^{1+2f(\theta)}} }[/math]

Аналогично целевой функции (1) значение связи между двумя кластерами делится на ожидаемое число связей между кластерами.

Вычислительное ядро алгоритма

Вычислительное ядро алгоритма состоит из двух частей:

1. Подсчет количества связей между точками;
2. Итеративное объединение кластеров.

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

Итеративная часть алгоритма может быть реализована с помощью структуры данных куча. На каждой итерации с каждым кластером [math]\displaystyle{ i }[/math] связана своя куча [math]\displaystyle{ q[i] }[/math], содержащая все связанный с ним кластеры, упорядоченная по убыванию функции качества. Также поддерживается глобальная куча [math]\displaystyle{ Q }[/math], содержащая все кластеры [math]\displaystyle{ j }[/math], упорядоченные по мере убывания функции качества [math]\displaystyle{ g(j, max(q[j])) }[/math], где [math]\displaystyle{ max(q[j]) }[/math] - это вершина кучи [math]\displaystyle{ q[j] }[/math]. Тогда на каждой итерации объединяются кластер [math]\displaystyle{ j }[/math] с вершины глобальной кучи и кластер [math]\displaystyle{ q[j] }[/math]. После объединения необходимо обновить все кучи, содержащие объединенные кластеры, обновить глобальную кучу [math]\displaystyle{ Q }[/math] и создать кучу для нового кластера. Итеративный процесс продолжается до тех пор, пока в глобальной куче [math]\displaystyle{ Q }[/math] не останется требуемое число кластеров.

Макроструктура алгоритма

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

Схема реализации последовательного алгоритма

Следующая схема последовательного алгоритма представлена в работе [1] (псевдокод):

procedure cluster(S, k)
begin
    link := compute_links(S)
    for s in S do
    begin
        q[s] := build_local_heap(link, s)
    end
    Q := build_global_heap(S, q)
    while size(Q) > k do
    begin
        u := extract_max(Q)
        v := max(q[u])
        delete(Q, v)
        w := merge(u, v)
        for each x in q[u] or q[v] do
        begin
            link[x, w] := link[x, u] + link[x, v]
            delete(q[x], u)
            delete(q[x], v)
            insert(q[x], w, g(x, w))
            insert(q[w], x, g(x, w))
            update(Q, x, q[x])
        end
        insert(Q, w, q[w])
        deallocate(q[u])
        deallocate(q[v])
    end
end

Последовательная сложность алгоритма

Для вычисления количества связей с помощью возведения матрицы в квадрат можно воспользоваться алгоритмом Копперсмита—Винограда [2], вычислительная сложность которого равна [math]\displaystyle{ O(n^{2,37}) }[/math] ([math]\displaystyle{ n }[/math] - количество кластеризуемых точек). Для случаев разреженной матрицы смежности в работе [1] описывается алгоритм, имеющий сложность [math]\displaystyle{ O(n^2 + n{m_m}{m_a}) }[/math], где [math]\displaystyle{ {m_a} }[/math] и [math]\displaystyle{ {m_m} }[/math] - это среднее и максимальное количество соседей точки соответственно.

Каждая куча может быть построена за время [math]\displaystyle{ O(n) }[/math] [3]. Внешний цикл итерационной части алгоритма исполняется [math]\displaystyle{ O(n) }[/math] раз. Так как на очередной итерации размер каждой локальной кучи в худшем случае может равняться [math]\displaystyle{ n }[/math], то новый кластер должен быть помещен в [math]\displaystyle{ n }[/math] локальных куч, что дает временную сложность вложенного цикла [math]\displaystyle{ O(n\log{n}) }[/math]. Таким образом, сложность внешнего цикла в худшем случае будет [math]\displaystyle{ O(n^2\log{n}) }[/math].

Итак, оценка времени выполнения алгоритма [math]\displaystyle{ C_t }[/math] задается следующей формулой:

[math]\displaystyle{ C_t = O(n^{2,37} + n^2\log{n}) }[/math]

Пространственную сложность подсчета количества связей можно оценить как [math]\displaystyle{ O(min(n^2, n{m_m}{m_a})) }[/math]. Так как при слиянии двух кластеров размер новой кучи не превосходит суммы размеров куч объединяемых кластеров, то пространственная сложность итерационной части алгоритма зависит только от размеров исходных куч и задается формулой:

[math]\displaystyle{ C_s = O(min(n^2, n{m_m}{m_a})) }[/math].

Информационный граф

Визуализация информационного графа для n = 100 исходных точек и k = 8 требуемых кластеров

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

• I - входные данные
• M - вычисление матрицы связей
• QS - созданиение локальной кучи
• QG - создание глобальной кучи
• HG - взятие элемента/элементов из кучи
• HD - удаление элемента из кучи
• HI - добавление элемента в кучу
• HU - обновление кучи
• HX - удаление куч
• CM - слияние кластеров
• A - арифметические операции
• K - полученный кластер

Ресурс параллелизма алгоритма

В описанном последовательном алгоритме основной ресурс параллелизма заключается в возможности параллельного вычисления количества связей (например, используя алгоритм параллельного вычисления произведения матриц), создания очередей, и параллельной обработки локальных очередей во внутреннем цикле итерационного этапа алгоритма. В таком случае при распараллеливании алгоритма на [math]\displaystyle{ p }[/math] вычислительных узлов верна следующая оценка времени работы алгоритма:

[math]\displaystyle{ C_{t_p} = {C_{L_p}} + O(\frac{n^2\log n}{p}) }[/math],

где [math]\displaystyle{ {C_{L_p}} }[/math] - оценка времени работы выбранного параллельного алгоритма подсчета числа связей.

Входные и выходные данные алгоритма

Входные данные:

1. Множество [math]\displaystyle{ n }[/math] точек, которые необходимо кластеризовать
2. Пороговое значение [math]\displaystyle{ \theta }[/math] из интервала [math]\displaystyle{ (0, 1) }[/math]
3. Требуемое число кластеров [math]\displaystyle{ k }[/math]

Также к входным данным относится информация о характеристиках точек (например в виде таблицы булевых значений принадлежности характеристик точкам).

Выходные данные:

1. [math]\displaystyle{ n }[/math] точек, распределенных по [math]\displaystyle{ k }[/math] кластерам

Свойства алгоритма

Рассматриваемый алгоритм является детерминированным.

Его вычислительная мощность равна

[math]\displaystyle{ \frac{O(n^{2,37} + n^2\log n)}{2n+2} }[/math]

Отношение последовательной и параллельной сложности:

[math]\displaystyle{ \frac{O(n^{2,37} + n^2\log n)}{{C_{L_p}} + O(\frac{n^2\log n}{p})} }[/math],

где [math]\displaystyle{ C_{L_p} }[/math] - сложность параллельного алгоритма подсчета числа связей (оценки приведены для использования алгоритма Копперсмита—Винограда для вычисления количества связей).

Одним из важных вопросов, влияющих на эффективность алгоритма, является выбор исходных точек. Подход к выбору размера исходного набора точек, достаточного для качественной кластеризации, описывается в работе [4]. Использование случайной стратегии выбора исходных точек помогает алгоритму кластеризации отбрасыванием выпадающих точек (то есть точек, слабо связанных с другими точками). Различные алгоритмы случайного выбора исходных точек рассматриваются в работе [5].

Обработка выпадающих точек в рассматриваемом алгоритме обеспечивается выбором порогового значения [math]\displaystyle{ \theta }[/math]. Также такие точки могут быть принудительно отброшены на последних итерациях алгоритма.

После завершения кластеризации выбранных исходных точек, оставшиеся точки распределяются по кластерам следующим образом. Из каждого кластера выбирается набор точек [math]\displaystyle{ {L_i} }[/math], и очередная точка [math]\displaystyle{ p }[/math] помещается в кластер [math]\displaystyle{ i }[/math], если [math]\displaystyle{ p }[/math] имеет наибольшее число соседей множеством [math]\displaystyle{ {L_i} }[/math].

Программная реализация алгоритма

Масштабируемость алгоритма и его реализации

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

Для оценки производительности была выбрана реализация из pyclustering, выполненная на C++. Для тестирования использовался суперкомпьютер «Ломоносов» [6] Суперкомпьютерного комплекса Московского университета.

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

Метод calculate_links в оригинальной версии отвечал за 20-25% общего времени выполнения. После распараллеливания на OpenMP увеличение количества узлов вело только к увеличению времени выполнения; оценки, полученные при распараллеливании этой функции далее не приводятся. Для наглядности исходный код метода приведен ниже:

size_t rock::calculate_links(const cluster & cluster1, const cluster & cluster2) const {
    size_t number_links = 0;
    #pragma omp parallel for reduction(+:number_links)
    for (size_t ii = 0; ii < cluster1.size(); ++ii) {
        for (size_t jj = 0; jj < cluster2.size(); ++jj) {
            std::size_t i = cluster1[ii];
            std::size_t j = cluster2[jj];
            number_links += (size_t) m_adjacency_matrix.has_connection(i, j);
        }
    }
}

Здесь для тестирования использовалась система с конфигурацией AMD FX-8120 4027 MHz, 24G DDR3 1866.

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

Метод create_adjacency_matrix также был распараллелен с помощью MPI. К качественному ускорению это не привело в первую очередь из-за того, что на использованных тестах среднее время выполнения этого метода занимало около 0.004% от общего времени выполнения. Распараллеленный код приведен ниже для справки.

void rock::create_adjacency_matrix(const dataset & p_data) {
    m_adjacency_matrix = adjacency_matrix(p_data.size());

    int procNum = 0;
    MPI_Comm_size(MPI_COMM_WORLD, &procNum);

    const size_t mSize = m_adjacency_matrix.size();
    const size_t q = mSize / procNum;
    const size_t blockSize = (mSize % procNum == 0) ? q : (q + 1);

    MPI_Request reqs[mSize];
    MPI_Request sendReqs[blockSize];
    MPI_Status stats[mSize];

    if (m_rank == 0) {
        for (size_t i = 0; i < blockSize; ++i) {
            reqs[i] = MPI_REQUEST_NULL;
        }
        for (size_t r = 1; r < procNum; ++r) {
            for (size_t i = r * blockSize; i < (r + 1) * blockSize && i < mSize; ++i) {
                MPI_Irecv(&(m_adjacency_matrix.raw_data()[i][0]), mSize, MPI_DOUBLE, r, 0, MPI_COMM_WORLD, reqs + i);
            }
        }
    }

    for (size_t i = m_rank * blockSize; i < (m_rank + 1) * blockSize && i < mSize; ++i) {
        for (size_t j = i + 1; j < mSize; j++) {
            double distance = euclidean_distance_sqrt(&p_data[i], &p_data[j]);
            if (distance < m_radius) {
                m_adjacency_matrix.set_connection(i, j);
            }
        }
    }

    if (m_rank == 0) {
        MPI_Waitall(mSize, reqs, stats);
    } else {
        size_t j = 0;
        for (size_t i = m_rank * blockSize; i < (m_rank + 1) * blockSize && i < mSize; ++i, ++j) {
            MPI_Isend(&(m_adjacency_matrix.raw_data()[i][0]), mSize, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD, sendReqs + j);
        }
        for (; j < blockSize; ++j) {
            sendReqs[j] = MPI_REQUEST_NULL;
        }
        MPI_Waitall(blockSize, sendReqs, stats);
    }
}

Тестирование было выполнено на кластере «Ломоносов» с узлами Intel Xeon E5-2697 v3 2.60GHz 64G DDR3. Результаты приведены ниже в форме графиков соответствия времени выполнения размеру задачи для последовательно выполняемой версии, а также зависимость времени выполнения тестов Chainlink (800 точек) и WingNut (1000 точек) от количества узлов.

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

Существующие реализации алгоритма

1. Библиотека annoviko/pyclustering включает реализации для Python и C++ под GPL 3
2. Пакет CBA из репозитория CRAN содержит реализацию на R под GPL 2
3. Библиотеки Weka и feed4weka содержат реализацию для Java под GPL 2

Литература

[1] Guha S., Rastogi R., Shim K. ROCK: A robust clustering algorithm for categorical attributes //Data Engineering, 1999. Proceedings., 15th International Conference on. – IEEE, 1999. – С. 512-521.

[2] Coppersmith D., Winograd S. Matrix multiplication via arithmetic progressions //Proceedings of the nineteenth annual ACM symposium on Theory of computing. – ACM, 1987. – С. 1-6.

[3] Thomas H. Cormen et al. Introduction to algorithms. – Cambridge : MIT press, 2001. – Т. 6.

[4] Guha S., Rastogi R., Shim K. Cure: an efficient clustering algorithm for large databases //Information Systems. – 2001. – Т. 26. – №. 1. – С. 35-58.

[5] Vitter J. S. Random sampling with a reservoir //ACM Transactions on Mathematical Software (TOMS). – 1985. – Т. 11. – №. 1. – С. 37-57.

[6] Воеводин Вл., Жуматий С., Соболев С., Антонов А., Брызгалов П., Никитенко Д., Стефанов К., Воеводин Вад. Практика суперкомпьютера «Ломоносов» // Открытые системы, 2012, N 7, С. 36-39.