Уровень алгоритма

Участник:Maria Zaitseva/PAM (Partitioning Around Medoids): различия между версиями

Материал из Алговики
Перейти к навигации Перейти к поиску
 
(не показаны 42 промежуточные версии 5 участников)
Строка 1: Строка 1:
 +
{{Assignment|ASA|Dexter}}
 +
 +
{{algorithm
 +
| name              = Partitioning Around Medoids
 +
| serial_complexity = <math>O(Tkn^2)</math>
 +
| input_data        = <math>\frac{n (n - 1)}{2}</math>
 +
| output_data      = <math>k</math>
 +
| pf_height        = <math>O(Tkn)</math>
 +
| pf_width          = <math>O(k(n-k))</math>
 +
}}
 +
 +
 +
 
Авторы: [[Участник: Maria Zaitseva|Зайцева М.Ф.]], [[Участник: kovalexal|Ковальчук А.А.]].
 
Авторы: [[Участник: Maria Zaitseva|Зайцева М.Ф.]], [[Участник: kovalexal|Ковальчук А.А.]].
  
Алгоритм кластеризации, основанный на разделении вокруг медоидов (Partitioning Around Medoids, PAM)
+
Авторы работали над каждым разделом статьи совместно (совместно искали и изучали материалы, формулировали описание и заполняли описание статьи на сайте). За содержание каждого из пунктов отвечают оба автора.   
 +
 
 +
'''PAM (Partitioning Around Medoids)''' – алгоритм разбиения множества объектов на [https://ru.wikipedia.org/wiki/Кластер кластеры], базирующийся на использовании наиболее представительных объектов каждого кластера ([https://ru.wikipedia.org/wiki/Медоид медоидов]) для их описания. Был разработан в 1987 году Леонардом Кауфманом и Питером Русивом <ref>Kaufman L., Rousseeuw P. J. Finding groups in data: an introduction to cluster analysis. – John Wiley & Sons, 2009. – Т. 344.</ref>. Показывает хорошие результаты разбиения на небольших наборах данных, однако имеет высокую вычислительную сложность, ввиду чего неэффективно работает на больших наборах данных.
  
= ЧАСТЬ. Свойства и структура алгоритмов =
+
= Свойства и структура алгоритма =
Свойства алгоритмов никак не зависят от вычислительных систем, и с этой точки зрения данная часть AlgoWiki имеет безусловную собственную ценность. Описание алгоритма делается один раз, после чего многократно используется для его реализации в различных программно-аппаратных средах. Несмотря на то, что в данной части мы рассматриваем лишь машинно-независимые свойства алгоритмов, соображения, важные на этапе реализации, или же ссылки на соответствующие пункты [[#ЧАСТЬ. Программная реализация алгоритмов|части II AlgoWiki]], здесь также вполне уместны.
 
  
 
== Общее описание алгоритма ==
 
== Общее описание алгоритма ==
В данном разделе представляется самый первый вариант описания тех задач (или классов задач), для решения которых предназначен алгоритм. В описании поясняются особенности как алгоритма, так и объектов, с которыми он работает. Если описание соответствует целому классу схожих по структуре алгоритмов, либо же посвящено описанию отдельного представителя некоторого класса, то это здесь указывается явно. Описываются базовые математические свойства и структура входных данных. При необходимости, в описании могут присутствовать формулы, а также даваться ссылки на описания других используемых алгоритмов. Данное описание должно быть достаточным для однозначного понимания сути решаемой задачи.
+
Данный алгоритм используется для решения задачи [https://ru.wikipedia.org/wiki/Кластерный_анализ кластеризации] множества объектов на заранее известное число кластеров <math>k</math>. Каждый кластер характеризуется наиболее представительным (центральным ) объектом – медоидом. Принадлежность объекта к соответствующему кластеру определяется индексом наиболее близкого к нему медоида.
 +
 
 +
На вход алгоритму может подаваться набор объектов с заданной функцией расчета дистанции между ними, либо матрица расстояний, характеризующая дистанции между каждой парой объектов, а также число кластеров, на которое необходимо разбить множество объектов. На выходе алгоритм выдает множество объектов, принятых в качестве медоидов.
 +
 
 +
Алгоритм состоит из двух последовательных этапов:
 +
 
 +
1. Шаг ''BUILD''. Построение начального набора <math>k</math> медоидов путем их последовательного выбора. Первый из них выбирается как объект, сумма дистанций которого до всех других объектов минимальна (центральный объект). Каждый следующий выбирается из оставшихся объектов как лучше всего минимизирующий сумму расстояний от объектов до ближайших к ним медоидов. Выбранный набор медоидов характеризует "конфигурацию" кластеров, для которой может быть вычислена ее стоимость.
 +
 
 +
2. Шаг ''SWAP''. На данном шаге производится попытка улучшить стоимость текущей конфигурации путем перебора всех возможных замен (медоид-немедоид). Среди рассматриваемых конфигураций ищется наиболее оптимальная. В случае, если данная конфигурация по стоимости лучше текущей, она принимается в качестве текущей и шаг повторяется с начала, иначе алгоритм завершается.
 +
 
 +
Преимущества данного алгоритма по сравнению с подобными ему:
 +
 
 +
* В отличие от более популярного алгоритма [https://ru.wikipedia.org/wiki/K-means k-means], PAM позволяет работать с произвольными объектами с заданной мерой рассчета расстояний (различий), не требует хранения всего множества объектов (достаточно знать лишь матрицу дистанций, так как медоиды выбираются из множества объектов), а также более устойчив к выбросам в данных (за счет рассмотрения медоидов вместо центроидов);
 +
* В отличие от традиционного [https://en.wikipedia.org/wiki/K-medoids k-medoids] на шаге ''BUILD'' начальные медоиды выбираются более оптимальным образом по сравнению со случайным выбором. Это позволяет уменьшить количество итераций шага ''SWAP'', тем самым сокращая время работы алгоритма.
 +
 
 +
Недостатком данного алгоритма является неэффективность при больших наборах данных ввиду необходимости полного перебора.
  
 
== Математическое описание алгоритма ==
 
== Математическое описание алгоритма ==
Приводится математическое описание решаемой задачи в виде совокупности формул и соотношений, как это принято в книгах и учебниках. По возможности, используются общепринятые обозначения и способы записи. Должны быть явно определены все использованные обозначения и описаны свойства входных данных. Представленное описание должно быть достаточным для однозначного понимания постановки решаемой задачи для человека, знающего математику.
+
''Входные данные'':
 +
 
 +
* <math>o = \{1,\dots,n\} \subset \mathbb{Z}</math> –– множество индексов объектов
 +
* <math>O = \{O_1,\dots,O_n\}</math> –– множество объектов
 +
* <math>k \leq n</math> –– количество кластеров для разбиения
 +
* <math>d:O\times O\rightarrow \mathbb{R}^{1}</math> –– метрика расчета расстояний между объектами, <math>d_{ij} = d(O_i,O_j)</math>
 +
* <math>D\in \mathbb{R}^{n\times n}</math> –– симметрическая матрица дистанций объектов, где <math>D=(d_{ij}), i \in o, j \in o</math>
 +
 
 +
''Выходные данные'':
 +
 
 +
* <math>m = \{m_1,\dots,m_k\}\subset o</math> –– множество индексов объектов, принятых в качестве медоидов; <math> M </math>--множество медоидов.
 +
 
 +
Принадлежность объекта к определенному кластеру определяется как индекс ближайшего к нему медоида:
 +
:<math>g_m(i)=\arg \min_{m_j\in m} d_{im_j}</math>.
 +
 
 +
Основная задача состоит в выборе такого множества индексов медоидов <math>m</math>, наилучшим образом минимизирующего функционал:
 +
:<math>F_o(m) = \sum_{j \in o \backslash m} d_{jg_m(j)}</math>.
 +
 
 +
Введем обозначения:
 +
 
 +
* <math>D_j = d_{jg_m(j)}</math> –– индекс первого по близости к j медоида
 +
* <math>E_j = d_{jg_{m \backslash g_m(j)}(j)}</math> –– индекс второго по близости к j медоида
 +
 
 +
Решение задачи состоит из следующих шагов:
 +
 
 +
1. Шаг ''BUILD''
 +
 
 +
  1.1. Выбор первого медоида
 +
:<math>m=\{ \arg \min_{i \in o} \sum_{j \in o} d_{ij} \}</math>
 +
 
 +
  1.2 Последовательный выбор оставшихся медоидов до тех пор, пока <math>\left\vert{m}\right\vert \neq k</math>
 +
:<math>m = m \cup \{ \arg \max_{i \in o \backslash m} \sum_{j \in o \backslash \{m \cup i\} } \max (D_j - d_{ji}, 0) \}</math>
 +
 
 +
2. Шаг ''SWAP''
 +
 
 +
  2.1. Производится перебор всех возможных пар медоид-немедоид <math>(i, h): i \in m, h \in o \backslash m</math>, для каждой из которых рассчитывается стоимость смены текущей конфигурации <math>T_m(i,h)</math>. При смене конфигурации множество медоидов изменяется <math>m = m \cup \{h\} \backslash \{i\}</math>.
 +
:<math>T_m(i, h) = \sum_{j \in o \backslash \{m \cup h\}} C_{jih}</math>, где <math>C_{jih}</math> –– вклад объекта <math>j</math> в перестановку медоида <math>i</math> и не медоида <math>h</math>.
 +
  При этом <math>C_{jih}</math> рассчитывается по-разному в зависимости от положения <math>j</math> относительно <math>i</math> и <math>h</math> (<math>j</math> дальше от <math>i</math> и от <math>h</math>, чем от своего ближайшего медоида; <math>i</math> оказался ближайшим медоидом <math>j</math>; <math>j</math> оказался дальше от <math>i</math>, но ближе к <math>h</math>, чем к своему медоиду).
 +
:<math>C_{jih} = \left\{\begin{matrix}
 +
0, & d_{ji} >  D_j, d_{jh} >  D_j \\
 +
d_{jh} - d_{ji}, & d_{ji} = D_j, d_{jh} < E_j \\
 +
E_j - D_j, & d_{ji} = D_j, d_{jh} \geqslant  E_j \\
 +
d_{jh} - D_j, & d_{ji} >  D_j, d_{jh} <  D_j
 +
\end{matrix}\right.</math>
 +
  Среди всех пар находится такая пара <math>(i_s, h_s)</math>, что <math>T_m(i_s, h_s)</math> минимально.
 +
 
 +
  2.2 В случае, если <math>T_m(i_s, h_s) \geq 0</math>, смена конфигурации медоидов не производится и решение задачи окончено. Иначе производится смена конфигурации и повторяется шаг 2.1.
  
 
== Вычислительное ядро алгоритма ==
 
== Вычислительное ядро алгоритма ==
В описываемом алгоритме выделяется и описывается [[глоссарий#Вычислительное ядро|''вычислительное ядро'']], т.е. та часть алгоритма, на которую приходится основное время работы алгоритма. Если в алгоритме несколько вычислительных ядер, то отдельно описывается каждое ядро. Описание может быть сделано в достаточно произвольной форме: словесной или с использованием языка математических формул. Вычислительное ядро может полностью совпадать с описываемым алгоритмом.
+
Вычислительное ядро алгоритма на этапе ''BUILD'' состоит последовательного выбора такого объекта из числа немедоидов, добавление которого к множеству медоидов максимизирует следующую функцию:
 +
 
 +
<math>F_{m} = \max_{i \in o \backslash m} \sum_{j \in o \backslash \{m \cup i\} } \max (D_j - d_{ji}, 0)</math>
 +
 
 +
С каждой следующей итерацией множество просматриваемых немедоидов уменьшается, таким образом, на <math>i</math>-ом шаге итерации необходимо рассчитывать сумму для <math>n - i</math> немедоидов.
 +
 
 +
Вычислительное ядро алгоритма на этапе ''SWAP'' состоит в расчете функции минимизации стоимости изменения конфигурации <math>T_m(i, h)</math> при замене медоида <math>i</math> на немедоид <math>h</math> в текущей конфигурации для всех возможных <math>(i, h): i \in m, h \in o \backslash m</math>:
 +
 
 +
<math>T_m(i, h) = \sum_{j \in o \backslash \{m \cup h\}} C_{jih}</math>.
  
 
== Макроструктура алгоритма ==
 
== Макроструктура алгоритма ==
Если алгоритм использует в качестве составных частей другие алгоритмы, то это указывается в данном разделе. Если в дальнейшем имеет смысл описывать алгоритм не в максимально детализированном виде (т.е. на уровне арифметических операций), а давать только его макроструктуру, то здесь описывается структура и состав макроопераций. Если в других разделах описания данного алгоритма в рамках AlgoWiki используются введенные здесь макрооперации, то здесь даются пояснения, необходимые для однозначной интерпретации материала. Типичные варианты макроопераций, часто встречающиеся на практике: нахождение суммы элементов вектора, скалярное произведение векторов, умножение  матрицы на вектор, решение системы линейных уравнений малого порядка, сортировка, вычисление значения функции в некоторой точке, поиск минимального значения в массиве, транспонирование матрицы, вычисление обратной матрицы и многие другие.
 
  
Описание макроструктуры очень полезно на практике. Параллельная структура алгоритмов может быть хорошо видна именно на макроуровне, в то время как максимально детальное отображение всех операций может сильно усложнить картину. Аналогичные аргументы касаются и многих вопросов реализации, и если для алгоритма эффективнее и/или технологичнее оставаться на макроуровне, оформив макровершину, например, в виде отдельной процедуры, то это и нужно отразить в данном разделе.  
+
1. Макрооперация "Выбор следующего медоида"
Выбор макроопераций не однозначен, причем, выделяя различные макрооперации, можно делать акценты на различных свойствах алгоритмов. С этой точки зрения, в описании одного алгоритма может быть представлено несколько вариантов его макроструктуры, дающих дополнительную информацию о его структуре. На практике, подобные альтернативные формы представления макроструктуры алгоритма могут оказаться исключительно полезными для его эффективной реализации на различных вычислительных платформах.
+
:<math>m = m \cup \{ \arg \max_{i \in o \backslash m} \sum_{j \in o \backslash \{m \cup i\} } \max (D_j - d_{ji}, 0) \}</math>
 +
 
 +
2. Макрооперация "Вычисление стоимости смены текущей конфигурации"
 +
:<math>T_m(i, h) = \sum_{j \in o \backslash \{m \cup h\}} C_{jih}</math>
 +
 
 +
Таким образом, алгоритм, в терминах макроопераций выглядит следующим образом:
 +
 
 +
1. Фаза ''BUILD''
 +
 
 +
На данном шаге осуществляется выбор первого медоида. Далее осуществляется макрооперация "Выбор следующего медоида" <math>(k - 1)</math> раз.
 +
 
 +
2. Фаза ''SWAP''
 +
 
 +
На данном шаге для всех возможных комбинаций медоид-немедоид осуществляется макрооперация "Вычисление стоимости смены текущей конфигурации". Данная операция производится <math>k(n-k)</math> раз для каждой итерации. Среди всех комбинаций выбирается та замена, стоимость которой максимальна.
 +
 
 +
3. ''Критерий останова''
 +
 
 +
:<math>T_m(i, h) < 0</math>
  
 
== Схема реализации последовательного алгоритма ==
 
== Схема реализации последовательного алгоритма ==
Здесь описываются все шаги, которые нужно выполнить при последовательной реализации данного алгоритма. В некотором смысле, данный раздел является избыточным, поскольку математическое описание уже содержит всю необходимую информацию. Однако он, несомненно, полезен: схема реализации алгоритма выписывается явно, помогая однозначной интерпретации приводимых далее оценок и свойств.
 
  
Описание может быть выполнено в виде блок-схемы, последовательности математических формул, обращений к описанию других алгоритмов, фрагмента кода на Фортране, Си или другом языке программирования, фрагмента кода на псевдокоде и т.п. Главное - это сделать схему реализации последовательного алгоритма полностью понятной. Совершенно не обязательно все шаги детализировать до элементарных операций, отдельные шаги могут соответствовать макрооперациям, отвечающим другим алгоритмам.
+
''Псевдокод алгоритма'':
  
Описание схемы реализации вполне может содержать и словесные пояснения, отражающие какие-либо тонкие нюансы самого алгоритма или его реализации. Уже в данном разделе можно сказать про возможный компромисс между объемом требуемой оперативной памяти и временем работы алгоритма, между используемыми структурами данных и степенью доступного параллелизма. В частности, часто возникает ситуация, когда можно ввести дополнительные временные массивы или же отказаться от использования специальных компактных схем хранения данных, увеличивая степень доступного параллелизма.
+
'''Вход''': Объем входных данных: <math>n\times n + 2</math>. Матрица дистанций D (<math>n\times n</math> вещественных чисел); количество кластеров <math>k \in \mathbb{Z}</math>; максимальное количество итераций алгоритма <math>maxIter \in \mathbb{Z}</math>
 +
'''Выход''': Объем выходных данных: <math>k</math>. Множество индексов медоидов <math>m = \{m_1, \dots, m_k\} \subset \mathbb{Z}</math>
 +
1 # Фаза BUILD
 +
2 Выбор стартового объекта в качестве начального медоида: m = {первый медоид}
 +
3 Пока |m| != k:
 +
4    Выбор следующего медоида i_m
 +
5    Добавление найденного медоида в множество: m = m + {i_m}
 +
6 # Фаза SWAP
 +
7 Пока (maxIter--):
 +
8    Лучшая стоимость замены: Tbest = MAXINT
 +
9    Кандидаты на замену: (i_s, h_s) = (-1, -1)
 +
10    Для всех i из множества медоидов:
 +
11        Для всех j из множества немедоидов:
 +
12            Расчет стоимости замены Tm(i,h)
 +
13            Если Tm(i, h) < Tbest:
 +
14                Tbest = Tm(i, j)
 +
15                Запомнить найденные кандидаты: (i_s, j_s) = (i, h)
 +
16    Если Tbest >= 0:
 +
17        Выход
 +
18    Иначе:
 +
19        m = m + {h} - {i}
 +
 
 +
На фазе ''BUILD'' осуществляется построение начального множества медоидов. Данный процесс производится итерационно до тех пор, пока не будет найдено необходимое число медоидов.
 +
На фазе ''SWAP'' производится попытка улучшить множество выбранных медоидов. Данный процесс также производится итерационно то тех пор, пока целевая функция не перестанет улучшаться. Также, для данного шага можно ввести ограничение по числу производимых алгоритмом итераций.
  
 
== Последовательная сложность алгоритма ==
 
== Последовательная сложность алгоритма ==
В данном разделе описания свойств алгоритма приводится оценка его [[глоссарий#Последовательная сложность|''последовательной сложности'']], т.е. числа операций, которые нужно выполнить при последовательном исполнении алгоритма (в соответствии с [[#Описание схемы реализации последовательного алгоритма|п.1.5]]). Для разных алгоритмов понятие операции, в терминах которой оценивается его сложность, может существенно различаться. Это могут быть операции для работы с вещественными числами, целыми числами, поразрядные операции, обращения в память, обновления элементов массива, элементарные функции, макрооперации и другие. В LU-разложении преобладают арифметические операции над вещественными числами, а для транспонирования матриц важны лишь обращения к памяти: это и должно найти отражение в описании.
 
  
Если выбор конкретного типа операций для оценки сложности алгоритма не очевиден, то нужно привести обоснование возможных вариантов. В некоторых случаях можно приводить оценку не всего алгоритма, а лишь его вычислительного ядра: в таком случае это нужно отметить, сославшись [[#Общее описание алгоритма|на п.1.1]].
+
Для осуществления шага ''BUILD'' потребуется:
  
Например, сложность алгоритма суммирования элементов вектора сдваиванием равна <math>n-1</math>. Сложность быстрого преобразования Фурье (базовый алгоритм Кули-Тьюки) для векторов с длиной, равной степени двойки – <math>n\log_2n</math> операций комплексного сложения и <math>(n\log_2n)/2</math> операций комплексного умножения. Сложность базового алгоритма разложения Холецкого (точечный вариант для плотной симметричной и положительно-определенной матрицы) это <math>n</math> вычислений квадратного корня, <math>n(n-1)/2</math> операций деления, по <math>(n^3-n)/6</math> операций умножения и сложения (вычитания).
+
* <math>n^2</math> операций сложения и <math>n</math> операций сравнения для выбора начального медоида;
 +
* <math>\sum_{l=1}^{k-1} (2(n-l) \cdot (n-l-1)+(n-l))=\frac{1}{6}(k-1)(4k^2 - 12kn + k + 6n(2n-1))</math> операций сложения, вычитания и сравнения для всех итераций выбора медоидов.
 +
 
 +
Таким образом, последовательная сложность шага ''BUILD'' оценивается как <math>O(kn^2)</math>.
 +
 
 +
Для осуществления каждого шага ''SWAP'' потребуется:
 +
 
 +
* <math>2k(n-k)(n-k-1)</math> операций вычитания и сложения;
 +
* <math>k(n-k)</math> операций сравнения;
 +
 
 +
Таким образом, последовательная сложность каждого шага ''SWAP'' оценивается как <math>O(kn^2)</math>.
 +
Заранее число итераций шага ''SWAP'' неизвестно, поэтому принимаем его равным <math>T</math>. Отсюда, общая последовательная сложность алгоритма PAM оценивается как <math>O(Tkn^2)</math>.
  
 
== Информационный граф ==
 
== Информационный граф ==
Это очень важный раздел описания. Именно здесь можно показать (увидеть) как устроена параллельная структура алгоритма, для чего приводится описание и изображение его информационного графа ([[глоссарий#Граф алгоритма|''графа алгоритма'']] [1]). Для рисунков с изображением графа будут составлены рекомендации по их формированию, чтобы все информационные графы, внесенные в энциклопедию, можно было бы воспринимать и интерпретировать одинаково. Дополнительно можно привести полное параметрическое  описание графа в терминах покрывающих функций [1].
 
  
Интересных вариантов для отражения информационной структуры алгоритмов много. Для каких-то алгоритмов нужно показать максимально подробную структуру, а иногда важнее макроструктура. Много информации несут разного рода проекции информационного графа, выделяя его регулярные составляющие и одновременно скрывая несущественные детали. Иногда оказывается полезным показать последовательность в изменении графа при изменении значений внешних переменных  (например, размеров матриц): мы часто ожидаем "подобное" изменение информационного графа, но это изменение не всегда очевидно на практике.
+
=== Информационный граф шага BUILD ===
 +
На рис.1 показана информационная структура шага BUILD алгоритма PAM на некотором шаге.
  
В целом, задача изображения графа алгоритма весьма нетривиальна. Начнем с того, что это потенциально бесконечный граф, число вершин и дуг которого определяется значениями внешних переменных, а они могут быть весьма и весьма велики. В такой ситуации, как правило, спасают упомянутые выше соображения подобия, делающие графы для разных значений внешних переменных "похожими": почти всегда достаточно привести лишь один граф небольшого размера, добавив, что графы для остальных значений будут устроены "точно также". На практике, увы, не всегда все так просто, и здесь нужно быть аккуратным.
+
[[file:Pam_build.png|thumb|center|553px|Рис.1. Информационная структура шага BUILD алгоритма PAM на шаге <math>t</math>]]
  
Далее, граф алгоритма - это потенциально многомерный объект. Наиболее естественная система координат для размещения вершин и дуг информационного графа опирается на структуру вложенности циклов в реализации алгоритма. Если глубина вложенности циклов не превышает трех, то и граф размещается в привычном трехмерном пространстве, однако для более сложных циклических конструкций с глубиной вложенности 4 и больше необходимы специальные методы представления и изображения графов.  
+
''Операции'':
 +
* <math>sum_{o,m_t}(i) = \sum_{j \in o \backslash \{m \cup i\}} \max (D_j, d_{ji})</math> –– оценка конфигурации при внесении <math>i</math> в множество индексов медоидов <math>m</math> с множеством объектов <math>o</math>;
 +
* <math>\arg \max</math> –– индекс объекта, дающего наибольшую оценку конфигурации при его внесении в множество индексов медоидов <math>m</math>;
 +
* <math>update(m)</math> –– операция добавления индекса объекта в множество индексов медоидов <math>m</math>.
  
В данном разделе AlgoWiki могут использоваться многие интересные возможности, которые еще подлежат обсуждению: возможность повернуть граф при его отображении на экране компьютера для выбора наиболее удобного угла обзора, разметка вершин по типу соответствующим им операций, отражение [[глоссарий#Ярусно-параллельная форма графа алгоритма|''ярусно-параллельной формы графа'']] и другие. Но в любом случае нужно не забывать главную задачу данного раздела - показать информационную структуру алгоритма так, чтобы стали понятны все его ключевые особенности, особенности параллельной структуры, особенности множеств дуг, участки регулярности и, напротив, участки с недерминированной структурой, зависящей от входных данных.  
+
=== Информационный граф шага SWAP ===
 +
На рис.2 показана информационная структура шага SWAP алгоритма PAM на некотором шаге.
  
На рис.1 показана информационная структура алгоритма умножения матриц, на рис.2 - информационная структура одного из вариантов алгоритма решения систем линейных алгебраических уравнений с блочно-двухдиагональной матрицей.
+
[[file:Pam_swap_new.png|thumb|center|553px|Рис.2. Информационная структура шага SWAP алгоритма PAM на шаге <math>t</math>]]
  
[[file:Fig1.svg|thumb|center|300px|Рис.1. Информационная структура алгоритма умножения матриц]]
+
''Операции'':
[[file:Fig2.svg|thumb|center|300px|Рис.2. Информационная структура одного из вариантов алгоритма решения систем линейных алгебраических уравнений с блочно-двухдиагональной матрицей]]
+
* <math>T_m(i, h) = \sum_{j \in o \backslash \{m \cup h\}} C_{jih} </math> –– стоимость смены конфигурации при замене немедоида <math>h</math> медоидом <math>i</math> в текущей конфигурации <math>m</math>;
 +
* <math>\arg \min</math> –– пара <math>(i, h)</math> (медоид-немедоид), замена которой в текущей конфигурации <math>m</math> дает наименьшую стоимость смены;
 +
* <math>swap</math> –– операция замены медоида <math>i</math> на немедоид <math>h</math> в текущей конфигурации <math>m</math>, если стоимость смены отрицательна.
  
 
== Ресурс параллелизма алгоритма ==
 
== Ресурс параллелизма алгоритма ==
Здесь приводится оценка [[глоссарий#Параллельная сложность|''параллельной сложности'']] алгоритма: числа шагов, за которое можно выполнить данный алгоритм в предположении доступности неограниченного числа необходимых процессоров (функциональных устройств, вычислительных узлов, ядер и т.п.). Параллельная сложность алгоритма понимается как высота канонической ярусно-параллельной формы [1]. Необходимо указать, в терминах каких операций дается оценка. Необходимо описать сбалансированность параллельных шагов по числу и типу операций, что определяется шириной ярусов канонической ярусно-параллельной формы и составом операций на ярусах.
+
Для осуществления шага ''BUILD'' потребуется:
  
Параллелизм в алгоритме часто имеет естественную иерархическую структуру. Этот факт очень полезен на практике, и его необходимо отразить в описании. Как правило, подобная иерархическая структура параллелизма хорошо отражается в последовательной реализации алгоритма через циклический профиль результирующей программы (конечно же, с учетом графа вызовов), поэтому циклический профиль ([[#Описание схемы реализации последовательного алгоритма|п.1.5]]) вполне  может быть использован и для отражения ресурса параллелизма.
+
* <math>n</math> параллельных операций сложения и <math>n</math> последовательных операций сравнения для выбора начального медоида;
 +
* <math>\sum_{l=1}^{k-1} 2(n-l) = (k-1)(2n-k)</math> параллельных операций сложения и вычитания для всех итераций выбора медоидов;
 +
* <math>\sum_{l=1}^{k-1} n-l = \frac{1}{2}(k-1)(2n-k)</math> последовательных операций сравнения для всех итераций выбора медоидов.
  
Для описания ресурса параллелизма алгоритма (ресурса параллелизма информационного графа) необходимо указать ключевые параллельные ветви в терминах [[глоссарий#Конечный параллелизм|''конечного'']] и [[глоссарий#Массовый параллелизм|''массового'']] параллелизма. Далеко не всегда ресурс параллелизма выражается просто, например, через [[глоссарий#Кооодинатный параллелизм|''координатный параллелизм'']] или, что то же самое, через независимость итераций некоторых циклов (да-да-да, циклы - это понятие, возникающее лишь на этапе реализации, но здесь все так связано… В данном случае, координатный параллелизм означает, что информационно независимые вершины лежат на гиперплоскостях, перпендикулярных одной из координатных осей). С этой точки зрения, не менее важен и ресурс [[глоссарий#Скошенный параллелизм|''скошенного параллелизма'']]. В отличие от координатного параллелизма, скошенный параллелизм намного сложнее использовать на практике, но знать о нем необходимо, поскольку иногда других вариантов и не остается: нужно оценить потенциал алгоритма, и лишь после этого, взвесив все альтернативы, принимать решение о конкретной параллельной реализации. Хорошей иллюстрацией может служить алгоритм, структура которого показана на рис.2: координатного параллелизма нет, но есть параллелизм скошенный, использование которого снижает сложность алгоритма с <math>n\times m</math> в последовательном случае до <math>(n+m-1)</math> в параллельном варианте.
+
Таким образом, параллельная сложность шага ''BUILD'' оценивается как <math>O(kn)</math>.
  
Рассмотрим алгоритмы, последовательная сложность которых уже оценивалась в [[#Последовательная сложность алгоритма|п.1.6]]. Параллельная сложность алгоритма суммирования элементов вектора сдваиванием равна <math>\log_2n</math>, причем число операций на каждом ярусе убывает с <math>n/2</math> до <math>1</math>. Параллельная сложность быстрого преобразования Фурье (базовый алгоритм Кули-Тьюки) для векторов с длиной, равной степени двойки - <math>\log_2n</math>. Параллельная сложность базового алгоритма разложения Холецкого (точечный вариант для плотной симметричной и положительно-определенной матрицы) это <math>n</math> шагов для вычислений квадратного корня, <math>(n-1)</math> шагов для операций деления и <math>(n-1)</math> шагов для операций умножения и сложения.
+
Для осуществления каждого шага ''SWAP'' потребуется:
 +
 
 +
* <math>k(n-k) </math> параллельных операций вычитания и сложения;
 +
* <math>k(n-k)</math> последовательных операций сравнения.
 +
 
 +
Таким образом, параллельная сложность каждого шага ''SWAP'' оценивается как <math>O(kn)</math>.
 +
Заранее число итераций шага ''SWAP'' неизвестно, поэтому принимаем его равным <math>T</math>. Отсюда, общая параллельная сложность алгоритма PAM оценивается как <math>O(Tkn)</math>.
  
 
== Входные и выходные данные алгоритма ==
 
== Входные и выходные данные алгоритма ==
В данном разделе необходимо описать объем, структуру, особенности и свойства входных и выходных данных алгоритма: векторы, матрицы, скаляры, множества, плотные или разреженные структуры данных, их объем. Полезны предположения относительно диапазона значений или структуры, например, диагональное преобладание в структуре входных матриц, соотношение между размером матриц по отдельным размерностям, большое число матриц очень малой размерности, близость каких-то значений к машинному нулю, характер разреженности матриц и другие.
+
 
 +
''Входные данные'':
 +
* Массив <math>V</math> из <math>\frac{n(n-1)}{2}</math> вещественных чисел, определяющий набор дистанций между всеми парами объектов (сжатая матрица дистанций). Симметрическая матрица дистанций  <math>D\in \mathbb{R}^{n\times n}</math>  восстанавливается из этого массива однозначным образом (для восстановления необходимо знать число объектов <math>n</math>):
 +
:<math>D(i, j)=\left\{\begin{matrix}
 +
V[ni - \frac{i(i + 1)}{2} + (j - i - 1)], & i < j\\
 +
V[nj - \frac{j(j + 1)}{2} + (i - j - 1)], & i > j \\
 +
0, & i = j
 +
\end{matrix}\right.</math>
 +
 
 +
* Целое число кластеров <math>k</math>, на которое необходимо разбить множество объектов.
 +
 
 +
''Выходные данные'':
 +
* Множество <math>m</math>, состоящее из <math>k</math> индексов объектов, принятых в качестве медоидов. Индекс кластера, которому принадлежит объект с индексом <math>i</math> можно восстановить с использованием функции <math>g_m(i)=\arg \min_{m_j\in m} d_{im_j}</math>.
  
 
== Свойства алгоритма ==
 
== Свойства алгоритма ==
Описываются прочие свойства алгоритма, на которые имеет смысл обратить внимание на этапе реализации. Как и ранее, никакой привязки к конкретной программно-аппаратной платформе не предполагается, однако вопросы реализации в проекте AlgoWiki всегда превалируют, и необходимость обсуждения каких-либо свойств алгоритмов определяется именно этим.
 
  
Весьма полезным является ''соотношение последовательной и параллельной сложности'' алгоритма. Оба понятия мы рассматривали ранее, но здесь делается акцент на том выигрыше, который теоретически может дать параллельная реализация алгоритма. Не менее важно описать и те сложности, которые могут возникнуть в процессе получения параллельной версии алгоритма.  
+
1. ''Соотношение последовательной и параллельной сложности'' в случае неограниченных ресурсов является ''линейным'' как отношение квадратичной сложности к линейной. 
 +
 
 +
2. ''Вычислительная мощность'' последовательного алгоритма, как отношение числа операций к суммарному объему входных и выходных данных: <math> O(Tk) </math>.  
 +
 
 +
3. Алгоритм является ''устойчивым'', так как все операции в алгоритме не подвержены накоплению ошибки.
 +
 
 +
4. Алгоритм является ''не детерминированным'', так как  не детерминирован выбор набора медоидов, на которых достигается наилучшая конфигурация.
  
[[глоссарий#Вычислительная мощность|''Вычислительная мощность'']] алгоритма равна отношению числа операций к суммарному объему входных и выходных данных. Она показывает, сколько операций приходится на единицу переданных данных. Несмотря на простоту данного понятия, это значение исключительно полезно на практике: чем выше  вычислительная мощность, тем меньше накладных расходов вызывает перемещение данных для их обработки, например, на сопроцессоре, ускорителе или другом узле кластера. Например, вычислительная мощность скалярного произведения двух векторов равна всего лишь <math>1</math>, а вычислительная мощность алгоритма умножения двух квадратных матриц равна <math>2n/3</math>.
+
5. Имеет высокую вычислительную сложность (квадратичную) ввиду необходимости полного перебора, поэтому неэффективно работает на больших наборах данных.
  
Вопрос первостепенной важности на последующем этапе реализации - это [[глоссарий#Устойчивость|''устойчивость'']] алгоритма. Все, что касается различных сторон этого понятия, в частности, оценки устойчивости, должно быть описано в данном разделе.
+
6. Позволяет работать с произвольными объектами с заданной мерой рассчета расстояний (различий), не требует хранения всего множества объектов.
  
''Сбалансированность'' вычислительного процесса можно рассматривать с разных сторон. Здесь и сбалансированность типов операций, в частности, арифметических операций между собой (сложение, умножение, деление) или же арифметических операций по отношению к операциям обращения к памяти (чтение/запись). Здесь и сбалансированность операций между параллельными ветвями алгоритма. С одной стороны, балансировка нагрузки является необходимым условием эффективной реализации алгоритма. Вместе с этим, это очень непростая задача, и в описании должно быть отмечено явно, насколько алгоритм обладает этой особенностью. Если обеспечение сбалансированности не очевидно, желательно описать возможные пути решения этой задачи.
+
7. Устойчив к выбросам данных (за счет рассмотрения медоидов).
  
На практике важна [[глоссарий#Детерминированность|''детерминированность алгоритмов'']], под которой будем понимать постоянство структуры вычислительного процесса. С этой точки зрения, классическое умножение плотных матриц является детерминированным алгоритмом, поскольку его структура при фиксированном размере матриц никак не зависит от элементов входных матриц. Умножение разреженных матриц, когда матрица хранятся в одном из специальных форматов, свойством детерминированности уже не обладает: его свойства, например, степень локальности данных зависит от структуры разреженности входных матриц. Итерационный алгоритм с выходом по точности также не является детерминированным: число итераций, а значит и число операций, меняется в зависимости от входных данных. В этом же ряду стоит использование датчиков случайных чисел, меняющих вычислительный процесс для различных запусков программы. Причина выделения свойства детерминированности понятна: работать с детерминированным алгоритмом проще, поскольку один раз найденная структура и будет определять качество его реализации. Если детерминированность нарушается, то это должно быть здесь описано вместе с описанием того, как недетерминированность влияет на структуру вычислительного процесса.
+
= Программная реализация алгоритма =
  
Серьезной причиной недетерминированности работы параллельных программ является изменение порядка выполнения ассоциативных операций. Типичный пример - это использование глобальных MPI-операций на множестве параллельных процессов, например, суммирование элементов распределенного массива. Система времени исполнения MPI сама выбирает порядок выполнения операций, предполагая выполнение свойства ассоциативности, из-за чего ошибки округления меняются от запуска программы к запуску, внося изменения в конечный результат ее работы. Это очень серьезная проблема, которая сегодня встречается часто на системах с массовым параллелизмом и определяет отсутствие повторяемости результатов работы параллельных программ. Данная особенность характерна для [[#ЧАСТЬ. Программная реализация алгоритмов|второй части AlgoWiki]], посвященной реализации алгоритмов, но вопрос очень важный, и соответствующие соображения, по возможности, должны быть отмечены и здесь.
+
== Особенности реализации последовательного алгоритма ==
  
Заметим, что, в некоторых случаях, недетерминированность в структуре алгоритмов можно "убрать" введением соответствующих макроопераций, после чего структура становится не только детерминированной, но и более понятной для восприятия. Подобное действие также следует отразить в данном разделе.
+
== Локальность данных и вычислений ==
  
[[глоссарий#Степень исхода|''Степень исхода вершины информационного графа'']] показывает, в скольких операциях ее результат будет использоваться в качестве аргумента. Если степень исхода вершины велика, то на этапе реализации алгоритма нужно позаботиться об эффективном доступе к результату ее работы. В этом смысле, особый интерес представляют рассылки данных, когда результат выполнения одной операции используется во многих других вершинах графа, причем число таких вершин растет с увеличением значения внешних переменных.
+
== Возможные способы и особенности параллельной реализации алгоритма ==
  
''"Длинные" дуги в информационном графе'' [1] говорят о потенциальных сложностях с размещением данных в иерархии памяти компьютера на этапе выполнения программы. С одной стороны, длина дуги зависит от выбора конкретной системы координат, в которой расположены вершины графа, а потому в другой системе координат они попросту могут исчезнуть (но не появится ли одновременно других длинных дуг?). А с другой стороны, вне зависимости от системы координат их присутствие может быть сигналом о необходимости длительного хранения данных на определенном уровне иерархии, что накладывает дополнительные ограничения на эффективность реализации алгоритма. Одной из причин возникновения длинных дуг являются рассылки скалярных величин по всем итерациям какого-либо цикла: в таком виде длинные дуги не вызывают каких-либо серьезных проблем на практике.
+
== Масштабируемость алгоритма и его реализации ==
  
Для проектирования специализированных процессоров или реализации алгоритма на ПЛИС представляют интерес ''компактные укладки информационного графа'' [1], которые также имеет смысл привести в данном разделе.
+
Для изучения масштабируемости алгоритма была выполнена его [https://gitlab.com/kovalexal/pam-lomonosov реализация] с использованием механизма MPI на языке программирования C++. Данная реализация использует библиотеку [http://www.boost.org Boost] (модуль MPI для удобного взаимодействия C++ и MPI и модуль uBLAS для работы с матрицами). Исследование масштабируемости описанной реализации производилось на суперкомпьютере [https://parallel.ru/cluster/lomonosov.html "Ломоносов"]. Для сборки программы использовался компилятор gcc версии 5.2.0 с опциями "-std=c++11 -O3" совместно с OpenMPI версии 1.8.4.
  
= ЧАСТЬ. Программная реализация алгоритма =
+
В качестве рассматриваемого при оценке набора данных была произведена их генерация на основе двумерного нормального распределения. Для этого сперва выбиралось некоторое количество точек, которые в дальнейшем выступали в качестве центроидов кластеров. Далее, для каждого полученного кластера осуществлялась генерация равного числа точек.
Вторая часть описания алгоритмов в рамках AlgoWiki рассматривает все составные части процесса их реализации. Рассматривается как последовательная реализация алгоритма, так и параллельная. Описывается взаимосвязь свойств программ, реализующих алгоритм, и особенностей архитектуры компьютера, на которой они выполняются. Исследуется работа с памятью, локальность данных и вычислений, описывается масштабируемость и эффективность параллельных программ, производительность компьютеров, достигаемая на данной программе. Обсуждаются особенности реализации для разных классов архитектур компьютеров, приводятся ссылки на реализации в существующих библиотеках.
 
  
== Особенности реализации последовательного алгоритма ==
 
Здесь описываются особенности и варианты реализации алгоритма в виде последовательной программы, которые влияют на [[глоссарий#Эффективность реализации|''эффективность ее выполнения'']]. В частности, в данном разделе имеет смысл ''сказать о существовании блочных вариантов реализации алгоритма'', дополнительно описав потенциальные преимущества или недостатки, сопровождающие такую реализацию. Важный вопрос - это ''возможные варианты организации работы с данными'', варианты структур данных, наборов временных массивов и другие подобные вопросы. Для различных вариантов реализации следует оценить доступный ресурс параллелизма и объем требуемой памяти.
 
  
Важным нюансом является ''описание необходимой разрядности выполнения операций алгоритма'' (точности). На практике часто нет никакой необходимости выполнять все арифметические операции над вещественными числами с двойной точностью, т.к. это не влияет ни на устойчивость алгоритма, ни на точность получаемого результата. В таком случае, если значительную часть операций можно выполнять над типом float, и лишь в некоторых фрагментах необходим переход к типу double, это обязательно нужно отметить. Это прямое указание не только на правильную реализацию с точки зрения устойчивости по отношению к ошибкам округления, но и на более эффективную.
+
[[file:points.png|thumb|center|553px|Рис.3. Пример сгенерированных данных – 20 кластеров, 5000 точек]]
  
Опираясь на информацию из [[#Описание ресурса параллелизма алгоритма|п.1.8]] (описание ресурса параллелизма алгоритма), при описании последовательной версии стоит сказать про возможности [[глоссарий#Эквивалентное преобразование|''эквивалентного преобразования программ'']], реализующих данных алгоритм. В дальнейшем, это даст возможность простого использования доступного параллелизма или же просто покажет, как использовать присущий алгоритму параллелизм на практике. Например, параллелизм на уровне итераций самого внутреннего цикла обычно используется для векторизации. Однако, в некоторых случаях этот параллелизм можно поднять "вверх" по структуре вложенности объемлющих циклов, что делает возможной и эффективную реализацию данного алгоритма на многоядерных SMP-компьютерах.
 
  
С этой же точки зрения, в данном разделе весьма полезны соображения по реализации алгоритма на различных параллельных вычислительных платформах. Высокопроизводительные кластеры, многоядерные узлы, возможности для векторизации или использования ускорителей - особенности этих архитектур не только опираются на разные свойства алгоритмов, но и по-разному должны быть выражены в программах, что также желательно описать в данном разделе.
+
Набор значений параметров запуска, используемых для оценке масштабируемости представлен ниже:
  
== [[Локальность данных и вычислений]] ==
+
* Число процессов: [1, 8, 16, 32, 48, 50, 51, 64, 99, 100, 101, 128, 160, 192, 200, 224, 256, 320, 384, 448, 512]
Вопросы локальности данных и вычислений не часто изучаются на практике, но именно локальность определяет эффективность выполнения программ на современных вычислительных платформах [2, 3]. В данном разделе приводятся оценки степени [[глоссарий#Локальность использования данных|''локальности данных'']] и [[глоссарий#Локальность вычислений|вычислений]] в программе, причем рассматривается как [[глоссарий#Временная локальность|''временна́я'']], так и [[глоссарий#Пространственная локальность|''пространственная'']] локальность. Отмечаются позитивные и негативные факты, связанные с локальностью, какие ситуации и при каких условиях могут возникать. Исследуется, как меняется локальность при переходе от последовательной реализации к параллельной. Выделяются ключевые шаблоны взаимодействия программы, реализующей описываемый алгоритм, с памятью. Отмечается возможная взаимосвязь между используемыми конструкциями языков программирования и степенью локальности, которыми обладают результирующие программы.
+
* Число кластеризуемых элементов: от 1000 до 5000 с шагом 250
 +
* Число кластеров при расчете: 20
  
Отдельно приводятся профили взаимодействия с памятью для вычислительных ядер и ключевых фрагментов. Если из-за большого числа обращений по общему профилю сложно понять реальную специфику взаимодействия программ с памятью, то проводится последовательная детализация и приводится серия профилей более мелкого масштаба.
 
  
На рис.3 и рис.4 показаны профили обращения в память для программ, реализующих разложение Холецкого и быстрое преобразование Фурье, по которым хорошо видна разница свойств локальности у данных алгоритмов.
+
Для оценки производительности алгоритма осуществлялся расчет числа операций с плавающей запятой в секунду (FLOPs или GFLOPs). Для этого при каждом запуске рассчитывалось время работы алгоритма, число итераций шага ''SWAP'', по данным значениям можно получить точное число операций сложения, вычитания и сравнения вещественных чисел (оно равно последовательной сложности алгоритма), поделив его на время работы можно получить оценку числа FLOP в секунду. График изменения производительности в зависимости от числа процессоров и числа кластеризуемых элементов представлен ниже. Интерактивный график доступен по [https://plot.ly/~kovalexal/7/ ссылке].
  
[[file:Cholesky_locality1.jpg|thumb|center|700px|Рис.3 Реализация метода Холецкого. Общий профиль обращений в память]]
+
[[file:PAM_perf.png|thumb|center|553px|Рис.4. Изменение производительности в зависимости от числа процессоров и числа кластеризуемых элементов]]
[[file:fft 1.PNG|thumb|center|700px|Рис.4 Нерекурсивная реализация БПФ для степеней двойки. Общий профиль обращений в память]]
 
  
== Возможные способы и особенности параллельной реализации алгоритма ==
+
Масштабируемость алгоритма была получена путем деления реальной производительности на пиковую произвдительность <math>s</math> задействованных ядер (пиковая производительность одного ядра "Ломоносова" составляет ''5.686'' GFLOPs). График изменения эффективности в зависимости от числа процессоров и числа кластеризуемых элементов представлен ниже. Интерактивный график доступен по [https://plot.ly/~kovalexal/9/ ссылке].
Раздел довольно обширный, в котором должны быть описаны основные факты и положения, формирующие параллельную программу. К их числу можно отнести:
 
* представленный иерархически ресурс параллелизма, опирающийся на структуру циклических конструкций и на граф вызовов программы;
 
* комбинацию (иерархию) массового параллелизма и параллелизма конечного;
 
* возможные способы распределения операций между процессами/нитями;
 
* возможные способы распределения данных;
 
* оценку количества операций, объёма и числа пересылок данных (как общего числа, так и в пересчёте на каждый параллельный процесс);
 
  
и другие.
+
[[file:PAM_eff.png|thumb|center|553px|Рис.5. Изменение эффективности в зависимости от числа процессоров и числа кластеризуемых элементов]]
  
В этом же разделе должны быть даны рекомендации или сделаны комментарии относительно реализации алгоритма с помощью различных технологий параллельного программирования: MPI, OpenMP, CUDA или использования директив векторизации.
 
  
== Масштабируемость алгоритма и его реализации ==
+
Из полученных значений эффективности можно оценить диапазон эффективности реализации алгоритма:
Задача данного раздела - показать пределы [[глоссарий#Масштабируемость|''масштабируемости'']] алгоритма на различных платформах. Очень важный раздел. Нужно выделить, описать и оценить влияние точек барьерной синхронизации, глобальных операций, операций сборки/разборки данных, привести оценки или провести исследование [[глоссарий#Сильная масштабируемость|''сильной'']] и [[глоссарий#Слабая масштабируемость|''слабой'']] масштабируемости алгоритма и его реализаций.
 
  
Масштабируемость алгоритма определяет свойства самого алгоритма безотносительно конкретных особенностей используемого компьютера. Она показывает, насколько параллельные свойства алгоритма позволяют использовать возможности растущего числа процессорных элементов. Масштабируемость параллельных программ определяется как относительно конкретного компьютера, так и относительно используемой технологии программирования, и в этом случае она показывает, насколько может вырасти реальная производительность данного компьютера на данной программе, записанной с помощью данной технологии программирования, при использовании бóльших вычислительных ресурсов (ядер, процессоров, вычислительных узлов).
+
* минимальная эффективность: 0.10353% достигается при использовании 512 процессов для кластеризации 1000 элементов на 20 кластеров;
 +
* максимальная эффективность: 0.47637% достигается при использовании 1 процесса для кластеризации 2000 элементов на 20 кластеров.
  
Ключевой момент данного раздела заключается в том, чтобы показать ''реальные параметры масштабируемости программы'' для данного алгоритма на различных вычислительных платформах в зависимости от числа процессоров и размера задачи  [4]. При этом важно подобрать такое соотношение между числом процессоров и размером задачи, чтобы отразить все характерные точки в поведении параллельной программы, в частности, достижение максимальной производительности, а также тонкие эффекты, возникающие, например, из-за блочной структуры алгоритма или иерархии памяти.
+
Построим оценки масштабируемости выбранной реализации алгоритма Partitioning Around Medoids согласно [[Scalability_methodology]]:
  
На рис.5. показана масштабируемость классического алгоритма умножения плотных матриц в зависимости от числа процессоров и размера задачи. На графике хорошо видны области с большей производительностью, отвечающие уровням кэш-памяти.
+
* По числу процессов: 6.07542E-04. При увеличении числа процессов эффективность реализации на рассмотренной области значений параметров запуска возрастает, но интенсивность данного роста мала. Для увеличения масштабируемости по числу процессов необходимо увеличить количество операций в секунду, осуществляемых каждым процессом (например, уменьшив время работы алгоритма). Для этого можно применить различные технологии распараллеливания и векторизации вычислений.
[[file:Масштабируемость перемножения матриц производительность.png|thumb|center|700px|Рис.5 Масштабируемость классического алгоритма умножения плотных матриц в зависимости от числа процессоров и размера задачи]]
+
* По размеру задачи: 1.91151E-04. При увеличении размера задачи эффективность возрастает, однако данный рост достаточно слаб. Это говорит нам о том, что эффективность алгоритма слабо зависит от размера входной задачи.  
 +
* По двум направлениям: 2.03868E-05. При рассмотрении увеличения как размера входа, так и числа процессов на всей рассмотренной области значений эффективность увеличивается, однако скорость увеличения эффективности небольшая. В совокупности с тем фактом, что разница между минимальной и максимальной эффективностью на рассмотренной области параметров небольшая, эффективность с увеличением масштабов возрастает, но медленно.
  
 
== Динамические характеристики и эффективность реализации алгоритма ==
 
== Динамические характеристики и эффективность реализации алгоритма ==
Это объемный раздел AlgoWiki, поскольку оценка эффективности реализации алгоритма требует комплексного подхода [5], предполагающего аккуратный анализ всех этапов от архитектуры компьютера до самого алгоритма. Основная задача данного раздела заключается в том, чтобы оценить степень эффективности параллельных программ, реализующих данный алгоритм на различных платформах, в зависимости от числа процессоров и размера задачи. Эффективность в данном разделе понимается широко: это и [[глоссарий#Эффективность распараллеливания|''эффективность распараллеливания'']] программы, это и [[глоссарий#Эффективность реализации|''эффективность реализации'']] программ по отношению к пиковым показателям работы вычислительных систем.
 
 
Помимо собственно показателей эффективности, нужно описать и все основные причины, из-за которых эффективность работы параллельной программы на конкретной вычислительной платформе не удается сделать выше. Это не самая простая задача, поскольку на данный момент нет общепринятой методики и соответствующего инструментария, с помощью которых подобный анализ можно было бы провести. Требуется оценить и описать эффективность работы с памятью (особенности профиля взаимодействия программы с памятью), эффективность использования заложенного в алгоритм ресурса параллелизма, эффективность использования коммуникационной сети (особенности коммуникационного профиля), эффективность операций ввода/вывода и т.п. Иногда достаточно интегральных характеристик по работе программы, в некоторых случаях полезно показать данные мониторинга нижнего уровня, например, по загрузке процессора, кэш-промахам, интенсивности использования сети Infiniband и т.п. Хорошее представление о работе параллельной MPI-программы дают данные трассировки, полученные, например, с помощью системы Scalasca.
 
  
 
== Выводы для классов архитектур ==
 
== Выводы для классов архитектур ==
В данный раздел должны быть включены рекомендации по реализации алгоритма для разных классов архитектур. Если архитектура какого-либо компьютера или платформы обладает специфическими особенностями, влияющими на эффективность реализации, то это здесь нужно отметить.
 
 
На практике это сделать можно по-разному: либо все свести в один текущий раздел, либо же соответствующие факты сразу включать в предшествующие разделы, где они обсуждаются и необходимы по смыслу. В некоторых случаях, имеет смысл делать отдельные варианты всей [[#ЧАСТЬ. Программная реализация алгоритмов|части II]] AlgoWiki применительно к отдельным классам архитектур, оставляя общей машинно-независимую [[#ЧАСТЬ. Свойства и структура алгоритмов|часть I]]. В любом случае, важно указать и позитивные, и негативные факты по отношению к конкретным классам. Можно говорить о возможных вариантах оптимизации или даже о "трюках" в написании программ, ориентированных на целевые классы архитектур.
 
  
 
== Существующие реализации алгоритма ==
 
== Существующие реализации алгоритма ==
Для многих пар алгоритм+компьютер уже созданы хорошие реализации, которыми можно и нужно пользоваться на практике. Данный раздел предназначен для того, чтобы дать ссылки на основные существующие последовательные и параллельные реализации алгоритма, доступные для использования уже сейчас. Указывается, является ли реализация коммерческой или свободной, под какой лицензией распространяется, приводится местоположение дистрибутива и имеющихся описаний. Если есть информация об особенностях, достоинствах и/или недостатках различных реализаций, то это также нужно здесь указать. Хорошими примерами реализации многих алгоритмов являются MKL, ScaLAPACK, PETSc, FFTW, ATLAS, Magma и другие подобные библиотеки.
 
  
= Литература =
+
1. Последовательные реализации
[1] Воеводин В.В., Воеводин Вл.В. Параллельные вычисления. - СПб.: БХВ-Петербург, 2002. - 608 с.
 
  
[2] Воеводин В.В., Воеводин Вад.В. Спасительная локальность суперкомпьютеров //Открытые системы. - 2013. - № 9. - С. 12-15.
+
* В [https://www.mathworks.com/help/stats/kmedoids.html?requestedDomain=www.mathworks.com#namevaluepairs MATLAB] ([https://ru.wikipedia.org/wiki/Коммерческое_программное_обеспечение Commercial software])
 +
* [http://elki.dbs.ifi.lmu.de/browser/elki/elki/src/main/java/de/lmu/ifi/dbs/elki/algorithm/clustering/kmeans/KMedoidsPAM.java ELKI] на языке JAVA (Свободное распространение ([https://ru.wikipedia.org/wiki/GNU_Affero_General_Public_License AGPL]))
 +
* В пакете [https://cran.r-project.org/package=cluster Cluster] на языке R (Свободное распространение ([https://ru.wikipedia.org/wiki/GNU_General_Public_License#GPL_v1 GPL]))
 +
* В пакете [https://cran.r-project.org/package=ClusterR ClusterR] на языке R (Свободное распространение ([https://ru.wikipedia.org/wiki/Лицензия_MIT MIT]))
  
[3] Воеводин Вад.В., Швец П. Метод покрытий для оценки локальности использования данных в программах // Вестник УГАТУ. — 2014. — Т. 18, № 1(62). — С. 224–229.
+
2. Параллельные реализации
  
[4] Антонов А.С., Теплов А.М. О практической сложности понятия масштабируемости параллельных программ// Высокопроизводительные параллельные вычисления на кластерных системах (HPC 2014): Материалы XIV Международной конференции -Пермь: Издательство ПНИПУ, 2014. С. 20-27.
+
* В пакете [https://cran.r-project.org/package=sprint SPRINT] на языке R (Свободное распространение ([https://ru.wikipedia.org/wiki/GNU_General_Public_License#GPL_v1 GPL]))
  
[5] Никитенко Д.А. Комплексный анализ производительности суперкомпьютерных систем, основанный на данных системного мониторинга // Вычислительные методы и программирование. 2014. 15. 85–97.
+
= Литература =
  
 
[[en:Description of algorithm properties and structure]]
 
[[en:Description of algorithm properties and structure]]

Текущая версия на 16:28, 19 декабря 2016

Symbol confirmed.svgЭта работа успешно выполнена
Преподавателю: в основное пространство, в подстраницу

Данное задание было проверено и зачтено.
Проверено Dexter и ASA.



Partitioning Around Medoids
Последовательный алгоритм
Последовательная сложность [math]O(Tkn^2)[/math]
Объём входных данных [math]\frac{n (n - 1)}{2}[/math]
Объём выходных данных [math]k[/math]
Параллельный алгоритм
Высота ярусно-параллельной формы [math]O(Tkn)[/math]
Ширина ярусно-параллельной формы [math]O(k(n-k))[/math]



Авторы: Зайцева М.Ф., Ковальчук А.А..

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

PAM (Partitioning Around Medoids) – алгоритм разбиения множества объектов на кластеры, базирующийся на использовании наиболее представительных объектов каждого кластера (медоидов) для их описания. Был разработан в 1987 году Леонардом Кауфманом и Питером Русивом [1]. Показывает хорошие результаты разбиения на небольших наборах данных, однако имеет высокую вычислительную сложность, ввиду чего неэффективно работает на больших наборах данных.

Содержание

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

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

Данный алгоритм используется для решения задачи кластеризации множества объектов на заранее известное число кластеров [math]k[/math]. Каждый кластер характеризуется наиболее представительным (центральным ) объектом – медоидом. Принадлежность объекта к соответствующему кластеру определяется индексом наиболее близкого к нему медоида.

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

Алгоритм состоит из двух последовательных этапов:

1. Шаг BUILD. Построение начального набора [math]k[/math] медоидов путем их последовательного выбора. Первый из них выбирается как объект, сумма дистанций которого до всех других объектов минимальна (центральный объект). Каждый следующий выбирается из оставшихся объектов как лучше всего минимизирующий сумму расстояний от объектов до ближайших к ним медоидов. Выбранный набор медоидов характеризует "конфигурацию" кластеров, для которой может быть вычислена ее стоимость.

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

Преимущества данного алгоритма по сравнению с подобными ему:

  • В отличие от более популярного алгоритма k-means, PAM позволяет работать с произвольными объектами с заданной мерой рассчета расстояний (различий), не требует хранения всего множества объектов (достаточно знать лишь матрицу дистанций, так как медоиды выбираются из множества объектов), а также более устойчив к выбросам в данных (за счет рассмотрения медоидов вместо центроидов);
  • В отличие от традиционного k-medoids на шаге BUILD начальные медоиды выбираются более оптимальным образом по сравнению со случайным выбором. Это позволяет уменьшить количество итераций шага SWAP, тем самым сокращая время работы алгоритма.

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

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

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

  • [math]o = \{1,\dots,n\} \subset \mathbb{Z}[/math] –– множество индексов объектов
  • [math]O = \{O_1,\dots,O_n\}[/math] –– множество объектов
  • [math]k \leq n[/math] –– количество кластеров для разбиения
  • [math]d:O\times O\rightarrow \mathbb{R}^{1}[/math] –– метрика расчета расстояний между объектами, [math]d_{ij} = d(O_i,O_j)[/math]
  • [math]D\in \mathbb{R}^{n\times n}[/math] –– симметрическая матрица дистанций объектов, где [math]D=(d_{ij}), i \in o, j \in o[/math]

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

  • [math]m = \{m_1,\dots,m_k\}\subset o[/math] –– множество индексов объектов, принятых в качестве медоидов; [math] M [/math]--множество медоидов.

Принадлежность объекта к определенному кластеру определяется как индекс ближайшего к нему медоида:

[math]g_m(i)=\arg \min_{m_j\in m} d_{im_j}[/math].

Основная задача состоит в выборе такого множества индексов медоидов [math]m[/math], наилучшим образом минимизирующего функционал:

[math]F_o(m) = \sum_{j \in o \backslash m} d_{jg_m(j)}[/math].

Введем обозначения:

  • [math]D_j = d_{jg_m(j)}[/math] –– индекс первого по близости к j медоида
  • [math]E_j = d_{jg_{m \backslash g_m(j)}(j)}[/math] –– индекс второго по близости к j медоида

Решение задачи состоит из следующих шагов:

1. Шаг BUILD

 1.1. Выбор первого медоида
[math]m=\{ \arg \min_{i \in o} \sum_{j \in o} d_{ij} \}[/math]
 1.2 Последовательный выбор оставшихся медоидов до тех пор, пока [math]\left\vert{m}\right\vert \neq k[/math]
[math]m = m \cup \{ \arg \max_{i \in o \backslash m} \sum_{j \in o \backslash \{m \cup i\} } \max (D_j - d_{ji}, 0) \}[/math]

2. Шаг SWAP

 2.1. Производится перебор всех возможных пар медоид-немедоид [math](i, h): i \in m, h \in o \backslash m[/math], для каждой из которых рассчитывается стоимость смены текущей конфигурации [math]T_m(i,h)[/math]. При смене конфигурации множество медоидов изменяется [math]m = m \cup \{h\} \backslash \{i\}[/math].
[math]T_m(i, h) = \sum_{j \in o \backslash \{m \cup h\}} C_{jih}[/math], где [math]C_{jih}[/math] –– вклад объекта [math]j[/math] в перестановку медоида [math]i[/math] и не медоида [math]h[/math].
 При этом [math]C_{jih}[/math] рассчитывается по-разному в зависимости от положения [math]j[/math] относительно [math]i[/math] и [math]h[/math] ([math]j[/math] дальше от [math]i[/math] и от [math]h[/math], чем от своего ближайшего медоида; [math]i[/math] оказался ближайшим медоидом [math]j[/math]; [math]j[/math] оказался дальше от [math]i[/math], но ближе к [math]h[/math], чем к своему медоиду).
[math]C_{jih} = \left\{\begin{matrix} 0, & d_{ji} \gt D_j, d_{jh} \gt D_j \\ d_{jh} - d_{ji}, & d_{ji} = D_j, d_{jh} \lt E_j \\ E_j - D_j, & d_{ji} = D_j, d_{jh} \geqslant E_j \\ d_{jh} - D_j, & d_{ji} \gt D_j, d_{jh} \lt D_j \end{matrix}\right.[/math]
 Среди всех пар находится такая пара [math](i_s, h_s)[/math], что [math]T_m(i_s, h_s)[/math] минимально.
 2.2 В случае, если [math]T_m(i_s, h_s) \geq 0[/math], смена конфигурации медоидов не производится и решение задачи окончено. Иначе производится смена конфигурации и повторяется шаг 2.1.

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

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

[math]F_{m} = \max_{i \in o \backslash m} \sum_{j \in o \backslash \{m \cup i\} } \max (D_j - d_{ji}, 0)[/math]

С каждой следующей итерацией множество просматриваемых немедоидов уменьшается, таким образом, на [math]i[/math]-ом шаге итерации необходимо рассчитывать сумму для [math]n - i[/math] немедоидов.

Вычислительное ядро алгоритма на этапе SWAP состоит в расчете функции минимизации стоимости изменения конфигурации [math]T_m(i, h)[/math] при замене медоида [math]i[/math] на немедоид [math]h[/math] в текущей конфигурации для всех возможных [math](i, h): i \in m, h \in o \backslash m[/math]:

[math]T_m(i, h) = \sum_{j \in o \backslash \{m \cup h\}} C_{jih}[/math].

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

1. Макрооперация "Выбор следующего медоида"

[math]m = m \cup \{ \arg \max_{i \in o \backslash m} \sum_{j \in o \backslash \{m \cup i\} } \max (D_j - d_{ji}, 0) \}[/math]

2. Макрооперация "Вычисление стоимости смены текущей конфигурации"

[math]T_m(i, h) = \sum_{j \in o \backslash \{m \cup h\}} C_{jih}[/math]

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

1. Фаза BUILD

На данном шаге осуществляется выбор первого медоида. Далее осуществляется макрооперация "Выбор следующего медоида" [math](k - 1)[/math] раз.

2. Фаза SWAP

На данном шаге для всех возможных комбинаций медоид-немедоид осуществляется макрооперация "Вычисление стоимости смены текущей конфигурации". Данная операция производится [math]k(n-k)[/math] раз для каждой итерации. Среди всех комбинаций выбирается та замена, стоимость которой максимальна.

3. Критерий останова

[math]T_m(i, h) \lt 0[/math]

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

Псевдокод алгоритма:

Вход: Объем входных данных: [math]n\times n + 2[/math]. Матрица дистанций D ([math]n\times n[/math] вещественных чисел); количество кластеров [math]k \in \mathbb{Z}[/math]; максимальное количество итераций алгоритма [math]maxIter \in \mathbb{Z}[/math]
Выход: Объем выходных данных: [math]k[/math]. Множество индексов медоидов [math]m = \{m_1, \dots, m_k\} \subset \mathbb{Z}[/math]
1 # Фаза BUILD
2 Выбор стартового объекта в качестве начального медоида: m = {первый медоид}
3 Пока |m| != k:
4     Выбор следующего медоида i_m
5     Добавление найденного медоида в множество: m = m + {i_m}
6 # Фаза SWAP
7 Пока (maxIter--):
8     Лучшая стоимость замены: Tbest = MAXINT
9     Кандидаты на замену: (i_s, h_s) = (-1, -1)
10     Для всех i из множества медоидов:
11         Для всех j из множества немедоидов:
12             Расчет стоимости замены Tm(i,h)
13             Если Tm(i, h) < Tbest:
14                 Tbest = Tm(i, j)
15                 Запомнить найденные кандидаты: (i_s, j_s) = (i, h)
16     Если Tbest >= 0:
17         Выход
18     Иначе:
19         m = m + {h} - {i}

На фазе BUILD осуществляется построение начального множества медоидов. Данный процесс производится итерационно до тех пор, пока не будет найдено необходимое число медоидов. На фазе SWAP производится попытка улучшить множество выбранных медоидов. Данный процесс также производится итерационно то тех пор, пока целевая функция не перестанет улучшаться. Также, для данного шага можно ввести ограничение по числу производимых алгоритмом итераций.

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

Для осуществления шага BUILD потребуется:

  • [math]n^2[/math] операций сложения и [math]n[/math] операций сравнения для выбора начального медоида;
  • [math]\sum_{l=1}^{k-1} (2(n-l) \cdot (n-l-1)+(n-l))=\frac{1}{6}(k-1)(4k^2 - 12kn + k + 6n(2n-1))[/math] операций сложения, вычитания и сравнения для всех итераций выбора медоидов.

Таким образом, последовательная сложность шага BUILD оценивается как [math]O(kn^2)[/math].

Для осуществления каждого шага SWAP потребуется:

  • [math]2k(n-k)(n-k-1)[/math] операций вычитания и сложения;
  • [math]k(n-k)[/math] операций сравнения;

Таким образом, последовательная сложность каждого шага SWAP оценивается как [math]O(kn^2)[/math]. Заранее число итераций шага SWAP неизвестно, поэтому принимаем его равным [math]T[/math]. Отсюда, общая последовательная сложность алгоритма PAM оценивается как [math]O(Tkn^2)[/math].

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

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

На рис.1 показана информационная структура шага BUILD алгоритма PAM на некотором шаге.

Рис.1. Информационная структура шага BUILD алгоритма PAM на шаге [math]t[/math]

Операции:

  • [math]sum_{o,m_t}(i) = \sum_{j \in o \backslash \{m \cup i\}} \max (D_j, d_{ji})[/math] –– оценка конфигурации при внесении [math]i[/math] в множество индексов медоидов [math]m[/math] с множеством объектов [math]o[/math];
  • [math]\arg \max[/math] –– индекс объекта, дающего наибольшую оценку конфигурации при его внесении в множество индексов медоидов [math]m[/math];
  • [math]update(m)[/math] –– операция добавления индекса объекта в множество индексов медоидов [math]m[/math].

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

На рис.2 показана информационная структура шага SWAP алгоритма PAM на некотором шаге.

Рис.2. Информационная структура шага SWAP алгоритма PAM на шаге [math]t[/math]

Операции:

  • [math]T_m(i, h) = \sum_{j \in o \backslash \{m \cup h\}} C_{jih} [/math] –– стоимость смены конфигурации при замене немедоида [math]h[/math] медоидом [math]i[/math] в текущей конфигурации [math]m[/math];
  • [math]\arg \min[/math] –– пара [math](i, h)[/math] (медоид-немедоид), замена которой в текущей конфигурации [math]m[/math] дает наименьшую стоимость смены;
  • [math]swap[/math] –– операция замены медоида [math]i[/math] на немедоид [math]h[/math] в текущей конфигурации [math]m[/math], если стоимость смены отрицательна.

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

Для осуществления шага BUILD потребуется:

  • [math]n[/math] параллельных операций сложения и [math]n[/math] последовательных операций сравнения для выбора начального медоида;
  • [math]\sum_{l=1}^{k-1} 2(n-l) = (k-1)(2n-k)[/math] параллельных операций сложения и вычитания для всех итераций выбора медоидов;
  • [math]\sum_{l=1}^{k-1} n-l = \frac{1}{2}(k-1)(2n-k)[/math] последовательных операций сравнения для всех итераций выбора медоидов.

Таким образом, параллельная сложность шага BUILD оценивается как [math]O(kn)[/math].

Для осуществления каждого шага SWAP потребуется:

  • [math]k(n-k) [/math] параллельных операций вычитания и сложения;
  • [math]k(n-k)[/math] последовательных операций сравнения.

Таким образом, параллельная сложность каждого шага SWAP оценивается как [math]O(kn)[/math]. Заранее число итераций шага SWAP неизвестно, поэтому принимаем его равным [math]T[/math]. Отсюда, общая параллельная сложность алгоритма PAM оценивается как [math]O(Tkn)[/math].

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

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

  • Массив [math]V[/math] из [math]\frac{n(n-1)}{2}[/math] вещественных чисел, определяющий набор дистанций между всеми парами объектов (сжатая матрица дистанций). Симметрическая матрица дистанций [math]D\in \mathbb{R}^{n\times n}[/math] восстанавливается из этого массива однозначным образом (для восстановления необходимо знать число объектов [math]n[/math]):
[math]D(i, j)=\left\{\begin{matrix} V[ni - \frac{i(i + 1)}{2} + (j - i - 1)], & i \lt j\\ V[nj - \frac{j(j + 1)}{2} + (i - j - 1)], & i \gt j \\ 0, & i = j \end{matrix}\right.[/math]
  • Целое число кластеров [math]k[/math], на которое необходимо разбить множество объектов.

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

  • Множество [math]m[/math], состоящее из [math]k[/math] индексов объектов, принятых в качестве медоидов. Индекс кластера, которому принадлежит объект с индексом [math]i[/math] можно восстановить с использованием функции [math]g_m(i)=\arg \min_{m_j\in m} d_{im_j}[/math].

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

1. Соотношение последовательной и параллельной сложности в случае неограниченных ресурсов является линейным как отношение квадратичной сложности к линейной.

2. Вычислительная мощность последовательного алгоритма, как отношение числа операций к суммарному объему входных и выходных данных: [math] O(Tk) [/math].

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

4. Алгоритм является не детерминированным, так как не детерминирован выбор набора медоидов, на которых достигается наилучшая конфигурация.

5. Имеет высокую вычислительную сложность (квадратичную) ввиду необходимости полного перебора, поэтому неэффективно работает на больших наборах данных.

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

7. Устойчив к выбросам данных (за счет рассмотрения медоидов).

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

2.1 Особенности реализации последовательного алгоритма

2.2 Локальность данных и вычислений

2.3 Возможные способы и особенности параллельной реализации алгоритма

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

Для изучения масштабируемости алгоритма была выполнена его реализация с использованием механизма MPI на языке программирования C++. Данная реализация использует библиотеку Boost (модуль MPI для удобного взаимодействия C++ и MPI и модуль uBLAS для работы с матрицами). Исследование масштабируемости описанной реализации производилось на суперкомпьютере "Ломоносов". Для сборки программы использовался компилятор gcc версии 5.2.0 с опциями "-std=c++11 -O3" совместно с OpenMPI версии 1.8.4.

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


Рис.3. Пример сгенерированных данных – 20 кластеров, 5000 точек


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

  • Число процессов: [1, 8, 16, 32, 48, 50, 51, 64, 99, 100, 101, 128, 160, 192, 200, 224, 256, 320, 384, 448, 512]
  • Число кластеризуемых элементов: от 1000 до 5000 с шагом 250
  • Число кластеров при расчете: 20


Для оценки производительности алгоритма осуществлялся расчет числа операций с плавающей запятой в секунду (FLOPs или GFLOPs). Для этого при каждом запуске рассчитывалось время работы алгоритма, число итераций шага SWAP, по данным значениям можно получить точное число операций сложения, вычитания и сравнения вещественных чисел (оно равно последовательной сложности алгоритма), поделив его на время работы можно получить оценку числа FLOP в секунду. График изменения производительности в зависимости от числа процессоров и числа кластеризуемых элементов представлен ниже. Интерактивный график доступен по ссылке.

Рис.4. Изменение производительности в зависимости от числа процессоров и числа кластеризуемых элементов

Масштабируемость алгоритма была получена путем деления реальной производительности на пиковую произвдительность [math]s[/math] задействованных ядер (пиковая производительность одного ядра "Ломоносова" составляет 5.686 GFLOPs). График изменения эффективности в зависимости от числа процессоров и числа кластеризуемых элементов представлен ниже. Интерактивный график доступен по ссылке.

Рис.5. Изменение эффективности в зависимости от числа процессоров и числа кластеризуемых элементов


Из полученных значений эффективности можно оценить диапазон эффективности реализации алгоритма:

  • минимальная эффективность: 0.10353% достигается при использовании 512 процессов для кластеризации 1000 элементов на 20 кластеров;
  • максимальная эффективность: 0.47637% достигается при использовании 1 процесса для кластеризации 2000 элементов на 20 кластеров.

Построим оценки масштабируемости выбранной реализации алгоритма Partitioning Around Medoids согласно Scalability_methodology:

  • По числу процессов: 6.07542E-04. При увеличении числа процессов эффективность реализации на рассмотренной области значений параметров запуска возрастает, но интенсивность данного роста мала. Для увеличения масштабируемости по числу процессов необходимо увеличить количество операций в секунду, осуществляемых каждым процессом (например, уменьшив время работы алгоритма). Для этого можно применить различные технологии распараллеливания и векторизации вычислений.
  • По размеру задачи: 1.91151E-04. При увеличении размера задачи эффективность возрастает, однако данный рост достаточно слаб. Это говорит нам о том, что эффективность алгоритма слабо зависит от размера входной задачи.
  • По двум направлениям: 2.03868E-05. При рассмотрении увеличения как размера входа, так и числа процессов на всей рассмотренной области значений эффективность увеличивается, однако скорость увеличения эффективности небольшая. В совокупности с тем фактом, что разница между минимальной и максимальной эффективностью на рассмотренной области параметров небольшая, эффективность с увеличением масштабов возрастает, но медленно.

2.5 Динамические характеристики и эффективность реализации алгоритма

2.6 Выводы для классов архитектур

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

1. Последовательные реализации

  • В MATLAB (Commercial software)
  • ELKI на языке JAVA (Свободное распространение (AGPL))
  • В пакете Cluster на языке R (Свободное распространение (GPL))
  • В пакете ClusterR на языке R (Свободное распространение (MIT))

2. Параллельные реализации

  • В пакете SPRINT на языке R (Свободное распространение (GPL))

3 Литература

  1. Kaufman L., Rousseeuw P. J. Finding groups in data: an introduction to cluster analysis. – John Wiley & Sons, 2009. – Т. 344.