Boruvka's algorithm

Алгоритм Борувки
Sequential algorithm
Serial complexity	[math]O(|E|ln(|V|))[/math]
Input data	[math]O(|V| + |E|)[/math]
Output data	[math]O(|V|)[/math]
Parallel algorithm
Parallel form height	[math]max O(ln(|V|)) [/math]
Parallel form width	[math]O(|E|)[/math]

1 Properties and structure of the algorithm

1.1 General description of the algorithm

Алгоритм Борувки предназначен для решения задачи о построении минимального остовного дерева во взвешенном неориентированном графе. Алгоритм хорошо параллелизуется и является основой для распределённого алгоритма GHS.

1.2 Mathematical description of the algorithm

Пусть задан связный неориентированный граф [math]G = (V, E)[/math] с весами рёбер [math]f(e)[/math]. Предполагается, что веса всех рёбер различны (если это не так, то можно упорядочить рёбра сначала по весу, а потом по номеру).

Алгоритм Борувки основан на следующих двух свойствах задачи:

• Минимальное ребро фрагмента. Пусть [math]F[/math] – фрагмент минимального остовного дерева и [math]e_F[/math] – ребро наименьшего веса, исходящее из [math]F[/math] (т.е. ровно один его конец является вершиной из [math]F[/math]). Если ребро [math]e_F[/math] единственно, то оно принадлежит минимальному остовному дереву.
• Схлопывание фрагментов. Пусть [math]F[/math] – фрагмент минимального остовного дерева графа [math]G[/math], а граф [math]G'[/math] получен из [math]G[/math] склеиванием вершин, принадлежащих [math]F[/math]. Тогда объединение [math]F[/math] и минимального остовного дерева графа [math]G'[/math] даёт минимальное остовное дерево исходного графа [math]G[/math].

В начале работы алгоритма каждая вершина графа [math]G[/math] является отдельным фрагментом. На очередном шаге у каждого фрагмента выбирается исходящее ребро минимального веса (если такое ребро существует). Выбранные рёбра добавляются в минимальное остовное дерево, а соответствующие фрагменты склеиваются.

1.3 Computational kernel of the algorithm

Основными операциями являются:

1. Поиск минимального по весу исходящего ребра в каждом фрагменте.
2. Объединение фрагментов.

1.4 Macro structure of the algorithm

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

Классический пример (из статьи Борувки) – спроектировать наиболее дешёвую электрическую сеть, зная стоимость устройства каждого участка электрической линии.

Пусть задан связный граф [math]G=(V,E)[/math] с вершинами [math]V = ( v_{1}, v_{2}, ..., v_{n} )[/math] и рёбрами [math]E = ( e_{1}, e_{2}, ..., e_{m} )[/math]. Каждому ребру [math]e \in E[/math] приписан вес [math]w(e)[/math].

Требуется построить дерево [math]T^* \subseteq E[/math], связывающее все вершины, и имеющее наименьший возможный вес среди всех таких деревьев:

[math] w(T^* )= \min_T( w(T)) [/math]

где вес множества рёбер есть сумма их весов:

[math]w(T)=\sum_{e \in T} (w(T))[/math]

Если граф [math]G[/math] не является связным, то дерева, связывающего все вершины, не существует.

В этом случае необходимо найти минимальной остовное дерево для каждой компоненты связности [math]G[/math]. Набор таких деревьев называется минимальным остовным лесом (сокращённо MSF – Minimum Spanning Forest).

1.4.1 Вспомогательный алгоритм: система непересекающихся множеств (Union-Find)

Во всех алгоритмах решения задачи требуется отслеживать, каким уже построенным фрагментам дерева принадлежат те или иные вершины графа. Для этого используется структура данных «система непересекающихся множеств» (Union-Find). Данная структура поддерживает две операции:

1. [math]FIND(v) = w[/math] – по вершине v возвращает вершину w – «корень» фрагмента, которому принадлежит вершина v. При этом гарантируется, что вершины u и v принадлежат одному и тому же фрагменту, тогда и только тогда, когда [math]FIND(u) = FIND(v)[/math].

2. [math]MERGE(u, v)[/math] – объединяет два фрагмента, которым принадлежат вершины [math]u[/math] и [math]v.[/math] (Если они уже лежат в одном фрагменте, то ничего не происходит.) При практической реализации удобно, чтобы данная операция возвращала значение истина, если объединение фрагментов имело место, и ложь в противном случае.

1.4.2 Последовательная версия

Классический последовательный алгоритм Union-Find описан в статье Тарьяна. Каждой вершине v приписывается указатель на вершину-родителя [math]parent(v)[/math].

1. Изначально [math]parent(v) := v[/math] для всех вершин.

2. [math]FIND(v)[/math] выполняется следующим образом: полагаем [math]u := v[/math], и далее следуем по указателям [math]u := parent(u)[/math] до тех пор, пока не станет [math]u = parent(u)[/math]. Это и будет результат операции. Дополнительно можно «схлопывать» дерево: присвоить всем посещённым вершинами: [math]parent(u_i) := u[/math], либо производить схлопывание по пути: [math]parent(u) := parent(parent(u)))[/math].

3. [math]MERGE(u, v)[/math] выполняется следующим образом: вначале находим корневые вершины [math]u := FIND(u), v := FIND(v)[/math]. Если [math]u = v[/math], то исходные вершины принадлежат одному фрагменту и объединения фрагментов не происходит. В противном случае полагаем одно из [math]parent(u) := v[/math] или [math]parent(v) := u[/math]. Дополнительно можно отслеживать количество вершин в каждом из фрагментов, чтобы меньший фрагмент подсоединять к большему, а не наоборот (оценки сложности доказываются именно при такой реализации, однако на практике алгоритм хорошо работает и без подсчёта количества вершин).

1.5 Implementation scheme of the serial algorithm

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

1. В самом начале каждая вершина является отдельным фрагментом.

2. На каждом шаге:

• Для каждого фрагмента определяется минимальное по весу исходящее ребро.

• Минимальные рёбра добавляются в минимальное остовное дерево, а соответствующие фрагменты объединяются.

3. Алгоритм останавливается, когда остаётся только один фрагмент, либо когда ни у одного из фрагментов нет исходящих рёбер.

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

Объединение фрагментов также может быть реализовано параллельно, с использованием описанной выше параллельной версии структуры Union-Find.

Аккуратный подсчёт количества активных фрагментов позволяет остановить алгоритм Борувки на один шаг раньше обычного:

1. В начале итерации счётчик активных фрагментов обнуляется.

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

3. На этапе объединения фрагментов счётчик уменьшается на единицу каждый раз, когда операция [math]MERGE(u, v)[/math] вернула значение истина.

Если в конце итерации счётчик равен 0 или 1, то вычисления останавливаются. Параллелизм возможен на этапе сортировки рёбер по весу, однако основной ход алгоритма является последовательным.

1.6 Serial complexity of the algorithm

Последовательная сложность алгоритма Борувки для графа с [math]|V|[/math] вершинами и [math]|E|[/math] рёбрами составляет [math]O(|E| \ln(|V|))[/math] операций.

1.7 Information graph

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

Нижний уровень параллелизма: поиск минимальных исходящих рёбер может выполняться независимо для каждого фрагмента, благодаря чему данную стадию вычислений можно эффективно параллелизовать (как на GPU, так и на CPU). Объединение фрагментов также может быть реализовано параллельно, с использованием описанной структуры Union-Find.

Верхний уровень параллелизма: построение отдельных минимальных основных деревьев для каждого из списков ребер может производиться параллельно. Например, список ребер может разбиваться на две части, одна из которых обрабатывается на GPU, а вторая параллельно на CPU.

Рисунок 1. Информационный граф нижнего уровня параллелизма

Рассмотрим информационные графы и подробное описание каждого из них. Так же можно считать, что на рисунке 1 представлен информационный граф классического алгоритма Борувки, а на рисунке 2 — алгоритма обработки графа.

Нижний уровень параллелизма на графе алгоритма (рисунок 1) расположен на уровнях {3, 4, 5}, соответствующим операциям параллельного поиска минимальных исходящих ребер, а так же уровнях {6, 7, 8}, соответствующим операциям параллельного объединения деревьев. Так же, различные операции копирования {1, 2, 8, 9} выполняются параллельно. После выполнения тела цикла, производится проверка {12} того, сколько деревьев осталось на текущем шаге, и если данное число не изменилось, то происходит выход из цикла, иначе аналогичная следующая итерация.

Верхний уровень параллелизма (рисунок 2), как уже говорилось, заключается в параллельном вычислении минимального основного дерева (compute mst) для различных частей графа. Перед этим производится процесс инициализации (init process), данные которого используют последующие параллельные compute mst. Затем, после параллельных вычислений mst, происходит вычисление итого основного дерево, после чего после чего полученный результат сохраняется (save results).

Рисунок 2. Информационный граф верхнего уровня параллелизма

1.8 Parallelization resource of the algorithm

Итак, алгоритм Борувки обладает двумя уровнями параллелизма.

На верхнем уровне минимальные остовные деревья могут искаться для отдельных частей списка ребер графа (параллельные compute_MST на рисунке 2). Однако, затем должно последовать финальное объединение полученных ребер и вычисление минимального остовного дерева для полученного графа, которое будет производиться последовательно.

Кроме того, вычисление каждого из минимальных остовных деревьев (параллельные compute_MST на рисунке 2) обладает внутренним ресурсом параллелизма, описанным далее. Операции инициализации и копирования данных ([1], [2], [9] на рисунке 1) могут производиться параллельно за [math]O(|V|)[/math] операций. Так же параллельно могут производиться операции поиска минимальных исходящих ребер ([3],[4],[5]), при том для каждой дуги и обратной к ней независимо, что даёт [math]2*O(|E|)[/math] параллельных операций. Помимо этого, операции объединения деревьев [6], [7], [8] могут так же производиться параллельно за [math]O(|V|)[/math] операций.

В результате, для классического алгоритма Борувки ширина ярусно-параллельной формы равна [math]O(|E|)[/math], а высота ЯПФ зависит от числа шагов алгоритма, и ограничена сверху [math]O(ln(|V|))[/math].

1.9 Input and output data of the algorithm

Входные данные: взвешенный граф [math](V, E, W)[/math] ([math]|V|[/math] вершин [math]v_i[/math] и [math]|E|[/math] рёбер [math]e_j = (v^{(1)}_{j}, v^{(2)}_{j})[/math] с весами [math]f_j[/math]).

Объём входных данных: [math]O(|V| + |E|)[/math].

Выходные данные: список рёбер минимального остовного дерева (для несвязного графа – список минимальных остовных деревьев для всех компонент связности).

Объём выходных данных: [math]O(|V|)[/math].

1.10 Properties of the algorithm

1. Алгоритм останавливается за конечное число шагов, поскольку на каждом шаге становится по крайней мере на один фрагмент меньше.
2. Более того, число фрагментов на каждом шаге уменьшается как минимум вдвое, так что общее число шагов составляет не более [math]\log_2 n[/math]. Отсюда следует и оценка сложности алгоритма.

2 Software implementation of the algorithm

2.2 Locality of data and computations

На этапе поиска минимального ребра происходят следующие обращения к памяти:

1. Чтение информации о рёбрах. Может производится последовательно.
2. Проверка принадлежности ребра одному и тому же фрагменту - два чтения массива [math]parent(u)[/math] с вероятным промахой по кэшу.
3. Чтение и обновление минимального веса ребра фрагмента. Данная информация может быть закеширована, особенно на поздних шагах, однако обновление необходимо производить атомарно, что требует инвализации кэша.

На этапе схлопывания фрагментов требуется атомарно обновить массив [math]parent(u)[/math] для каждого добавляемого в MST ребра. В зависимости от реализации параллельной структуры Union-Find корни фрагментов могут находиться ближе к началу массива, что позволяет закешировать эту наиболее часто читаемую область. Требование атомарности, однако, ограничивает эффект от такого кэширования.

2.2.1 Locality of implementation

2.2.1.1 Structure of memory access and a qualitative estimation of locality

Рисунок 3. Алгоритм Борувки. Общий профиль обращений в память

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

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

Начнем с изучения верхней части профиля (фрагмент 1 на рис. 3), которая показана на рис. 4. На каждом этапе большую часть обращений занимает последовательный перебор всех элементов данного фрагмента (выделен на рис. 4 желтым). Остальные обращения на разных этапах устроены по-разному. Если на первом этапе эти обращения разбросаны достаточно далеко друг от друга, что приводит к низкой пространственной и временной локальности, то на последнем этапе почти все обращения (не считая последовательного перебора) выполняются к одному и тому же элементу, что, естественно, характеризуется очень высокой локальностью. Подобное строение всего фрагмента приводит, скорее всего, к средним значениям и по пространственной, и по временной локальности.

Рисунок 4. Профиль обращений, фрагмент 1

Далее перейдем к изучению фрагмента 2 (рис. 5). Здесь можно увидеть, что строение каждого из 4 этапов отличается достаточно сильно. Как и в случае с фрагментом 1, каждый следующий этап обладает более высокой локальностью, однако здесь это заметно сильнее. При этом отметим, что данный фрагмент задействует всего около 60 элементов, а обращений к ним выполняется достаточно много, так что локальность в данном случае будет высока.

Рисунок 5. Профиль обращений, фрагмент 2

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

Рисунок 6. Профиль обращений, фрагмент 3

Отдельное рассмотрение фрагмента 4 (рис. 7) позволяет увидеть, что локальность здесь определяется 4 последовательными переборами всех элементов данного фрагмента. Эти переборы обладают стандартной структурой – шаг по памяти 1, только 1 обращение к каждому элементу; небольшое искривление данных переборов вызвано нерегулярной активностью в других фрагментах, которая приводит к искажению визуального представления профиля. Подобный набор обращений обладает высокой пространственной, но низкой временной локальностью.

Рисунок 7. Профиль обращений, фрагмент 4

Таким образом, фрагменты 2 и 3 характеризуются высокой локальностью, другие 2 фрагмента – средней локальностью. А поскольку большая часть обращений приходится именно на фрагменты 2 и 3, можно предположить, что общая локальность должна быть достаточно высока.

2.2.1.2 Quantitative estimation of locality

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

На рисунке 8 приведены значения daps для реализаций распространенных алгоритмов, отсортированные по возрастанию (чем больше daps, тем в общем случае выше производительность). Можно увидеть, что производительность работы с памятью в данном случае достаточно неплоха – значение daps сравнимо, например, со значением для реализации метода Холецкого. Однако это значение заметно ниже самых производительных реализаций алгоритмов (например, теста Linpack), что в целом неудивительно в случае графовых алгоритмов, традиционно неэффективно работающих с памятью.

Рисунок 8. Сравнение значений оценки daps

2.3 Possible methods and considerations for parallel implementation of the algorithm

Программа, реализующая алгоритм Борувки, состоит из двух частей:

1. части, отвечающей за общую координацию вычислений

2. части, отвечающей за параллельные вычисления на многоядерных CPU или GPU.

Описанный выше последовательный алгоритм не может применяться в параллельной программе: в реализации [math]MERGE[/math] результаты операций [math]FIND(u)[/math] и [math]FIND(v)[/math] могут постоянно меняться, что приведёт к race condition. Параллельный вариант алгоритма описан в статье.

1. Каждой вершине v соответствует запись [math]A[v] = { parent, rank }[/math]. Изначально [math]A[v] := { v, 0 }[/math].

2. Вспомогательная операция [math]UPDATE(v, rank_v, u, rank_u)[/math]:

old := A[v]

if old.parent != v or old.rank != rank_v then return false

new := { u, rank_u }

return CAS(A[v], old, new)

3. Операция [math]FIND(v)[/math]:

   while v != A[v].parent do
       u := A[v].parent
       CAS(A[v].parent, u, A[u].parent)
       v := A[u].parent
   return v

4. Операция UNION(u, v):

   while true do
       (u, v) := (FIND(u), FIND(v))
       if u = v then return false
       (rank_u, rank_v) := (A[u].rank, A[v].rank)
       if (A[u].rank, u) > (A[v].rank, v) then
           swap((u, rank_u), (v, rank_v))
       if UPDATE(u, rank_u, v, rank_u) then
           if rank_u = rank_v then
               UPDATE(v, rank_v, v, rank_v + 1)
           return true

Для описанной версии алгоритма гарантируется свойство wait-free. На практике может использоваться упрощённая версия без подсчёта рангов, обладающая более слабым свойством lock-free, но в ряде случаев выигрывающая по скорости.

2.4 Scalability of the algorithm and its implementations

2.4.1 Scalability of the algorithm

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

1. пропускная способность памяти при чтении данных графа
2. соперничество потоков при выполнении атомарных операций с памятью
3. барьерная синхронизация после каждого подшага алгоритма.

2.4.2 Scalability of of the algorithm implementation

Проведём исследование масштабируемости параллельной реализации алгоритма Борувки согласно методике. Исследование проводилось на суперкомпьютере "Ломоносов-2 Суперкомпьютерного комплекса Московского университета.

Набор и границы значений изменяемых параметров запуска реализации алгоритма:

• число процессоров [1 : 28] с шагом 1;
• размер графа [2^20 : 2^27].

Проведем отдельные исследования сильной масштабируемости и масштабируемости вширь реализации алгоритма Борувки.

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

Рисунок 3. Параллельная реализация алгоритма Борувки масштабируемость CPU версии: производительность в зависимости от числа запущенных CPU-потоков.

Рисунок 4. Параллельная реализация алгоритма Борувки масштабируемость различных версий реализации алгоритма: производительность в зависимости от размера графа

2.5 Dynamic characteristics and efficiency of the algorithm implementation

Для проведения экспериментов использовалась реализация алгоритма Борувки, реализованная для CPU. Все результаты получены на суперкомпьютере «Ломоносов-2». Использовались процессоры Intel Xeon E5-2697v3, задача решалась для графа большого размера на одном узле. На рисунках показана эффективность реализации алгоритма Борувки, запуск проводился на 1 узле для графа 2^27, выполнялась 1 итерация.

Рисунок 9. График загрузки CPU при выполнении алгоритма Борувки

На графике загрузки процессора видно, что почти все время работы программы не загружены и средний уровень загрузки составляет около 10%. Это нормальная картина для программ, запущенных c использованием одного ядра в реализации.

Рисунок 10. График числа процессов, ожидающих вхождения в стадию счета (Loadavg), при работе алгоритма Борувки

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

Рисунок 11. График кэш-промахов L1 в секунду при работе алгоритма Борувки

На графике кэш-промахов первого уровня видно, что число промахов очень высокое и находится на уровне 40 млн/сек . Интересен факт увеличения числа промахов к концу итерации до уровня в 70 млн/сек.

Рисунок 12. График кэш-промахов L1 в секунду при работе алгоритма Борувки

На графике кэш-промахов второго уровня видно, что число промахов достаточно тоже высокое и находится на уровне 30 млн/сек . На графике промахов второго уровня факт увеличения числа промахов к концу итерации проявляется более явно и увеличивается до значения в 50млн/сек.

Рисунок 13. График кэш-промахов L3 в секунду при работе алгоритма Борувки

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

Рисунок 13. График скорости передачи по сети Infiniband в байт/сек при работе алгоритма Борувки

На графике скорости передачи данных по сети Infiniband наблюдается достаточно высокая интенсивность использования сети на кпервом этапе. Это объясняется программной логикой, которая предполагает чтение графа из файла с диска, коммуникации с которым происходят на Ломоносов-2 через выделенную для этих задач сеть Infiniband.

Рисунок 14. График скорости передачи по сети Infiniband в пакетах/сек при работе алгоритма Борувки

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

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

2.7 Existing implementations of the algorithm

• C++, MPI: Parallel Boost Graph Library; функции dense_boruvka_minimum_spanning_tree, boruvka_then_merge, boruvka_mixed_merge сочетают алгоритм Борувки и алгоритм Крускала.