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

Участник:Sergey Lavrushkin/EM-алгоритм кластеризации: различия между версиями

Материал из Алговики
Перейти к навигации Перейти к поиску
 
(не показано 80 промежуточных версий 3 участников)
Строка 1: Строка 1:
 +
{{Assignment|Algoman|Dexter}}
 +
{{algorithm
 +
| name              = EM-алгоритм кластеризации (1 итерация)
 +
| serial_complexity = <math>O(nkq^2)</math>
 +
| pf_height        = <math>O(\log{nk})</math>
 +
| pf_width          = <math>O(nk)</math>
 +
| input_data        = <math>nq + 1</math>
 +
| output_data      = <math>k(n + q^2 + q + 1)</math>
 +
}}
 
Автор статьи: [[Участник: Sergey Lavrushkin|Сергей Лаврушкин]] (группа 620)
 
Автор статьи: [[Участник: Sergey Lavrushkin|Сергей Лаврушкин]] (группа 620)
  
Строка 6: Строка 15:
  
 
==== EM-алгоритм ====
 
==== EM-алгоритм ====
EM-алгоритм (англ. expectation-maximization) — алгоритм, используемый в математической статистике для нахождения оценок максимального правдоподобия параметров вероятностных моделей, в случае, когда модель зависит от некоторых скрытых переменных. Каждая итерация алгоритма состоит из двух шагов. На E-шаге (expectation) вычисляется ожидаемое значение функции правдоподобия, при этом скрытые переменные рассматриваются как наблюдаемые. На M-шаге (maximization) вычисляется оценка максимального правдоподобия, таким образом увеличивается ожидаемое правдоподобие, вычисляемое на E-шаге. Затем это значение используется для E-шага на следующей итерации. Алгоритм выполняется до сходимости.<ref>ЕМ-алгоритм, его модификации и обобщения [http://www.machinelearning.ru/wiki/index.php?title=EM-%D0%B0%D0%BB%D0%B3%D0%BE%D1%80%D0%B8%D1%82%D0%BC http://www.machinelearning.ru/]</ref>
+
EM-алгоритм (англ. expectation-maximization) — алгоритм, используемый в математической статистике для нахождения оценок максимального правдоподобия параметров вероятностных моделей, в случае, когда модель зависит от некоторых скрытых переменных. Каждая итерация алгоритма состоит из двух шагов. На E-шаге (expectation) вычисляется ожидаемое значение функции правдоподобия, при этом скрытые переменные рассматриваются как наблюдаемые. На M-шаге (maximization) вычисляется оценка максимального правдоподобия, таким образом увеличивается ожидаемое правдоподобие, вычисляемое на E-шаге. Затем это значение используется для E-шага на следующей итерации. Алгоритм выполняется до сходимости.<ref>ЕМ-алгоритм, его модификации и обобщения [http://www.machinelearning.ru/wiki/index.php?title=EM-%D0%B0%D0%BB%D0%B3%D0%BE%D1%80%D0%B8%D1%82%D0%BC http://www.machinelearning.ru/wiki/index.php?title=EM-алгоритм]</ref>
  
 
Как правило, ЕМ-алгоритм применяется для решения задач двух типов:
 
Как правило, ЕМ-алгоритм применяется для решения задач двух типов:
Строка 16: Строка 25:
 
==== Кластеризация с помощью EM-алгоритма ====
 
==== Кластеризация с помощью EM-алгоритма ====
  
Кластеризация — задача разбиения заданной выборки объектов (ситуаций) на непересекающиеся подмножества, называемые кластерами, так, чтобы каждый кластер состоял из схожих объектов, а объекты разных кластеров существенно отличались <ref>Кластеризация [http://www.machinelearning.ru/wiki/index.php?title=%D0%9A%D0%BB%D0%B0%D1%81%D1%82%D0%B5%D1%80%D0%B8%D0%B7%D0%B0%D1%86%D0%B8%D1%8F www.machinelearning.ru]</ref>. В основе идеи применения EM-алгоритма для кластеризации лежит предположение, что исследуемое множество данных может быть смоделировано с помощью линейной комбинации многомерных нормальных распределений, а целью является оценка параметров распределения, которые максимизируют логарифмическую функцию правдоподобия, используемую в качестве меры качества модели. К оцениваемым параметрам относятся математические ожидания и ковариационные матрицы. Предполагается, что любое наблюдение принадлежит ко всем кластерам, но с разной вероятностью. Задача заключается в "подгонке" распределений смеси к данным, а затем в определении вероятностей принадлежности наблюдения к каждому кластеру. Наблюдение должно быть отнесено к тому кластеру, для которого данная вероятность выше.
+
Кластеризация — задача разбиения заданной выборки объектов (ситуаций) на непересекающиеся подмножества, называемые кластерами, так, чтобы каждый кластер состоял из схожих объектов, а объекты разных кластеров существенно отличались <ref>Кластеризация [http://www.machinelearning.ru/wiki/index.php?title=%D0%9A%D0%BB%D0%B0%D1%81%D1%82%D0%B5%D1%80%D0%B8%D0%B7%D0%B0%D1%86%D0%B8%D1%8F http://www.machinelearning.ru/wiki/index.php?title=Кластеризация]</ref>. В основе идеи применения EM-алгоритма для кластеризации лежит предположение, что исследуемое множество данных может быть смоделировано с помощью линейной комбинации многомерных нормальных распределений, а целью является оценка параметров распределения, которые максимизируют логарифмическую функцию правдоподобия, используемую в качестве меры качества модели. К оцениваемым параметрам относятся математические ожидания и ковариационные матрицы. Предполагается, что любое наблюдение принадлежит ко всем кластерам, но с разной вероятностью. Задача заключается в "подгонке" распределений смеси к данным, а затем в определении вероятностей принадлежности наблюдения к каждому кластеру. Наблюдение должно быть отнесено к тому кластеру, для которого данная вероятность выше.
  
Алгоритм EM основан на вычислении расстояний. Он может рассматриваться как обобщение кластеризации на основе анализа смеси вероятностных распределений. В процессе работы алгоритма происходит итеративное улучшение решения, а остановка осуществляется в момент, когда достигается требуемый уровень точности модели. Мерой в данном случае является монотонно увеличивающаяся статистическая величина, называемая логарифмическим правдоподобием.<ref>EМ — масштабируемый алгоритм кластеризации [https://basegroup.ru/community/articles/em basegroup.ru]</ref>
+
Алгоритм EM основан на вычислении расстояний. Он может рассматриваться как обобщение кластеризации на основе анализа смеси вероятностных распределений. В процессе работы алгоритма происходит итеративное улучшение решения, а остановка осуществляется в момент, когда достигается требуемый уровень точности модели. Мерой в данном случае является монотонно увеличивающаяся статистическая величина, называемая логарифмическим правдоподобием.<ref>EМ — масштабируемый алгоритм кластеризации [https://basegroup.ru/community/articles/em https://basegroup.ru/community/articles/em]</ref>
  
 
=== Математическое описание алгоритма ===
 
=== Математическое описание алгоритма ===
 
Исходные данные:  
 
Исходные данные:  
 
* <math>k</math> — число кластеров,
 
* <math>k</math> — число кластеров,
* <math>X =\{x_{1}, x_{2}, ..., x_{n}\}</math> — множество из <math>n</math> наблюдений <math>q</math>-мерного пространства.
+
* <math>X =\{x_{1}, x_{2}, ..., x_{n}\}</math> — множество из <math>n</math> наблюдений размерности <math>q</math>: <math>x_{i} = (x_{i1}, ..., x_{iq})</math>.
 
Вычисляемые данные:
 
Вычисляемые данные:
 
* <math>(w, \theta) = (w_1, ..., w_k; \mu_1, ..., \mu_k; \Sigma_1, ..., \Sigma_k)</math> — параметры смеси гауссовых распределений,
 
* <math>(w, \theta) = (w_1, ..., w_k; \mu_1, ..., \mu_k; \Sigma_1, ..., \Sigma_k)</math> — параметры смеси гауссовых распределений,
Строка 43: Строка 52:
 
</math>
 
</math>
 
где <math>\theta_j</math> — параметры распределения <math>p_j(x)</math>.
 
где <math>\theta_j</math> — параметры распределения <math>p_j(x)</math>.
В случае смести гауссовых распределений предполагается, что все компоненты имеют многомерное нормальное распределение. То есть в случае <math>q</math>-мерного пространства признаков <math>\phi(\theta_j; x)</math> имеют следующий вид:
+
В случае смеси гауссовых распределений в <math>q</math>-мерном пространстве признаков <math>\phi(\theta_j; x)</math> имеют следующий вид:
 
:<math>
 
:<math>
 
\begin{align}
 
\begin{align}
& \phi(\theta_j; x) = \frac{1}{\Big(2\pi\Big)^{\frac{q}{2}}\sqrt{|\Sigma_j|}}\exp{\bigg\{-\frac{\delta^2}{2}\bigg\} },\\
+
\phi(\theta_j; x) = \frac{1}{\Big(2\pi\Big)^{\frac{q}{2}}\sqrt{|\Sigma_j|}}\exp{\bigg\{-\frac{\delta^2}{2}\bigg\} },  
& \delta^2 = \Big( x - \mu_j \Big)^T \Sigma_j^{-1} \Big( x - \mu_j \Big),
 
 
\end{align}
 
\end{align}
 
</math>
 
</math>
Строка 53: Строка 61:
 
* <math>\Sigma_j</math> — ковариационная матрица размером <math>q \times q</math>,
 
* <math>\Sigma_j</math> — ковариационная матрица размером <math>q \times q</math>,
 
* <math>\mu_j</math> — <math>q</math>-мерный вектор математических ожиданий,
 
* <math>\mu_j</math> — <math>q</math>-мерный вектор математических ожиданий,
* <math>\delta^2</math> — квадратичное расстояние Махаланобиса.
+
* <math>\delta^2 = \Big( x - \mu_j \Big)^T \Sigma_j^{-1} \Big( x - \mu_j \Big)</math> — квадратичное расстояние Махаланобиса.
 
В случае смеси гауссовых распределений получаем, что <math>\theta_j = (\mu_j, \Sigma_j)</math>. То есть для решения задачи разделения смеси гауссовых распределений необходимо оценить вектор параметров <math>(w, \theta) = (w_1, ..., w_k; \mu_1, ..., \mu_k; \Sigma_1, ..., \Sigma_k)</math>. Согласно принципу максимизации правдоподобия:
 
В случае смеси гауссовых распределений получаем, что <math>\theta_j = (\mu_j, \Sigma_j)</math>. То есть для решения задачи разделения смеси гауссовых распределений необходимо оценить вектор параметров <math>(w, \theta) = (w_1, ..., w_k; \mu_1, ..., \mu_k; \Sigma_1, ..., \Sigma_k)</math>. Согласно принципу максимизации правдоподобия:
 
:<math>
 
:<math>
Строка 60: Строка 68:
 
\end{align}
 
\end{align}
 
</math>
 
</math>
Таким образом, имеет место задача максимизации суммы логарифмов сумм, решение которой представляет большую трудность. В таком случае полезным оказывается итеративный метод решения — EM-алгоритм.<ref>Поиск структуры в данных — Lecture 9 — Expectation Maximization (EM-алгоритм) [https://www.coursera.org/learn/unsupervised-learning/lecture/2jhbI/expectation-maximization-em-alghoritm www.coursera.org]</ref>
+
Таким образом, имеет место задача максимизации суммы логарифмов сумм, решение которой представляет большую трудность. В таком случае полезным оказывается итеративный метод решения — EM-алгоритм.<ref>Поиск структуры в данных — Lecture 9 — Expectation Maximization (EM-алгоритм) [https://www.coursera.org/learn/unsupervised-learning/lecture/2jhbI/expectation-maximization-em-alghoritm https://www.coursera.org/learn/unsupervised-learning/lecture/2jhbI/expectation-maximization-em-alghoritm]</ref>
  
 
==== EM-алгоритм ====
 
==== EM-алгоритм ====
Оптимальные параметры задачи разделения смеси гауссовых распределений отыскиваются последовательно с помощью итерационного EM-алгоритма. Основная идея – вводится вспомогательный вектор скрытых переменных. Это позволяет свести сложную оптимизационную задачу к последовательности итераций по пересчету коэффициентов (скрытых переменных по текущему приближению вектора параметров - E-шаг) и максимизации правдоподобия (с целью найти следующее приближение вектора - М-шаг).<ref>EM алгоритм (пример) [http://www.machinelearning.ru/wiki/index.php?title=EM_%D0%B0%D0%BB%D0%B3%D0%BE%D1%80%D0%B8%D1%82%D0%BC_(%D0%BF%D1%80%D0%B8%D0%BC%D0%B5%D1%80) www.machinelearning.ru]</ref>
+
Оптимальные параметры задачи разделения смеси гауссовых распределений отыскиваются последовательно с помощью итерационного EM-алгоритма. Основная идея – вводится вспомогательный вектор скрытых переменных. Это позволяет свести сложную оптимизационную задачу к последовательности итераций по пересчету коэффициентов (скрытых переменных по текущему приближению вектора параметров - E-шаг) и максимизации правдоподобия (с целью найти следующее приближение вектора - М-шаг).<ref>EM алгоритм (пример) [http://www.machinelearning.ru/wiki/index.php?title=EM_%D0%B0%D0%BB%D0%B3%D0%BE%D1%80%D0%B8%D1%82%D0%BC_(%D0%BF%D1%80%D0%B8%D0%BC%D0%B5%D1%80) http://www.machinelearning.ru/wiki/index.php?title=EM_алгоритм_(пример)]</ref>
  
 
EM-алгоритм заключается в следующем:
 
EM-алгоритм заключается в следующем:
 
<ol>
 
<ol>
<li> В начале работы алгоритма задаются параметры начального приближения. Наиболее общим способом инициализации является присвоение элементам матрицы математических ожиданий случайных значений <math> \mu_j \leftarrow Random </math>, начальные ковариационные матрицы определяются как единчиные <math>\Sigma_j \leftarrow I</math>, веса кластеров задаются одинаковыми <math>w_i \leftarrow \frac{1}{k}</math>. Также в качестве начальных параметров можно использовать результат работы алгоритма K-means. (Данная эвристика применяется, так как K-means требуется намного меньше итераций до достижения стабилизации, в то время как каждый шаг EM требует больших вычислительных затрат).</li>
+
<li> В начале работы алгоритма задаются параметры начального приближения. Наиболее общим способом инициализации является присвоение элементам матрицы математических ожиданий случайных значений <math> \mu_j \leftarrow Random </math>, начальные ковариационные матрицы определяются как единичные <math>\Sigma_j \leftarrow I</math>, веса кластеров задаются одинаковыми <math>w_i \leftarrow \frac{1}{k}</math>. Также в качестве начальных параметров можно использовать результат работы алгоритма K-means. (Данная эвристика применяется, так как K-means требуется намного меньше итераций до достижения стабилизации, в то время как каждый шаг EM требует больших вычислительных затрат).</li>
 
<li> Далее итеративно выполняется следующая пара процедур:
 
<li> Далее итеративно выполняется следующая пара процедур:
 
<ul>
 
<ul>
<li>E-шаг: используя текущее значение вектора параметров <math>(w, \theta) = (w_1, ..., w_k; \mu_1, ..., \mu_k; \Sigma_1, ..., \Sigma_k)</math>, вычисляем значение вектора скрытых переменных <math>y</math>:
+
<li>E-шаг: используя текущее значение вектора параметров <math>(w, \theta) = (w_1, ..., w_k; \mu_1, ..., \mu_k; \Sigma_1, ..., \Sigma_k)</math>, вычисляются значения скрытых переменных <math>y_{ij}</math>:
 
:<math>
 
:<math>
 
\begin{align}
 
\begin{align}
Строка 76: Строка 84:
 
\end{align}
 
\end{align}
 
</math></li>
 
</math></li>
<li>М-шаг: переоценка вектора параметров, используя текущее значение вектора скрытых переменных:
+
<li>М-шаг: переоценка вектора параметров, используя текущие значения скрытых переменных:
 
:<math>
 
:<math>
 
\begin{align}
 
\begin{align}
Строка 112: Строка 120:
 
         <math>sump_i = 0</math>
 
         <math>sump_i = 0</math>
 
         '''for''' <math>j = \overline{1, k}</math>:
 
         '''for''' <math>j = \overline{1, k}</math>:
             <math>\delta_{ij} = \big( x_i - \mu_j \big)^T \Sigma_j^{-1} \big( x_i - \mu_j \big)</math>
+
             <math>\delta_{ij}^2 = \big( x_i - \mu_j \big)^T \Sigma_j^{-1} \big( x_i - \mu_j \big)</math>
             <math>g_{ij} = \frac{w_j}{{2 \pi}^\frac{q}{2}|\Sigma_j|}\exp{\big\{-\frac{\delta^2}{2}\big\} }</math>
+
             <math>g_{ij} = \frac{w_j}{{(2 \pi)}^\frac{q}{2}\sqrt{|\Sigma_j|}}\exp{\big\{-\frac{\delta_{ij}^2}{2}\big\} }</math>
 
             <math>sump_i \mathrel{+}= g_{ij}</math>
 
             <math>sump_i \mathrel{+}= g_{ij}</math>
 
         <math>y_i = \frac{g_i}{sump_i}</math>
 
         <math>y_i = \frac{g_i}{sump_i}</math>
Строка 179: Строка 187:
 
Граф алгоритма состоит из четырех групп вершин:
 
Граф алгоритма состоит из четырех групп вершин:
 
*'''Первая''' группа вершин соответствует вычислению обратных матриц и определителей ковариационных матриц. Число вершин: <math>k</math>. Аргументом этой операции является ковариационная матрица <math>\Sigma_j</math>. Результат срабатывания операции является промежуточным данным алгоритма.
 
*'''Первая''' группа вершин соответствует вычислению обратных матриц и определителей ковариационных матриц. Число вершин: <math>k</math>. Аргументом этой операции является ковариационная матрица <math>\Sigma_j</math>. Результат срабатывания операции является промежуточным данным алгоритма.
*'''Вторая''' группа вершин соответствует вычислению <math>g_{ij} = \frac{w_j}{{2 \pi}^\frac{q}{2}|\Sigma_j|}\exp{\big\{-\frac{1}{2}  \big( x_i - \mu_j \big)^T \Sigma_j^{-1} \big( x_i - \mu_j \big)  \big\}}</math>. Число вершин: <math>nk</math>. Аргументами этой операции являются вес кластера <math>w_j</math>, математическое ожидание <math>\mu_j</math>, ковариационная матрица <math>\Sigma_j</math>, наблюдение элемент выборки <math>x_i</math>. Результат срабатывания операции является промежуточным данным алгоритма.
+
*'''Вторая''' группа вершин соответствует вычислению <math>g_{ij} = \frac{w_j}{{2 \pi}^\frac{q}{2}|\Sigma_j|}\exp{\big\{-\frac{1}{2}  \big( x_i - \mu_j \big)^T \Sigma_j^{-1} \big( x_i - \mu_j \big)  \big\}}</math>. Число вершин: <math>nk</math>. Аргументами этой операции являются вес кластера <math>w_j</math>, математическое ожидание <math>\mu_j</math>, ковариационная матрица <math>\Sigma_j</math>, элемент выборки <math>x_i</math>. Результат срабатывания операции является промежуточным данным алгоритма.
 
*'''Третья''' группа вершин соответствует вычислению нормировочных коэффициентов <math>\sum_{l=1}^{k}g_{ij}</math>. Число вершин: <math>n</math>. Аргументами этой операции являются значения <math>g_{ij}, \, j = \overline{1, k}</math>. Результат срабатывания операции является промежуточным данным алгоритма.
 
*'''Третья''' группа вершин соответствует вычислению нормировочных коэффициентов <math>\sum_{l=1}^{k}g_{ij}</math>. Число вершин: <math>n</math>. Аргументами этой операции являются значения <math>g_{ij}, \, j = \overline{1, k}</math>. Результат срабатывания операции является промежуточным данным алгоритма.
 
*'''Четвертая''' группа вершин соответствует вычислению скрытых переменных <math>y_{ij}</math>. Число вершин: <math>nk</math>. Аргументами этой операции являются значение <math>g_{ij}</math> и нормировочный коэффициент. Результат срабатывания операции является выходным данным алгоритма.
 
*'''Четвертая''' группа вершин соответствует вычислению скрытых переменных <math>y_{ij}</math>. Число вершин: <math>nk</math>. Аргументами этой операции являются значение <math>g_{ij}</math> и нормировочный коэффициент. Результат срабатывания операции является выходным данным алгоритма.
Строка 208: Строка 216:
  
 
=== Ресурс параллелизма алгоритма ===
 
=== Ресурс параллелизма алгоритма ===
 +
 +
По графам из пункта 1.7 видно, что скрытые переменные <math>y_{ij}</math> в E-шаге, <math>y_{ij}x_i</math> в вычислении математических ожиданий <math>\mu_{j}</math> и <math>y_{ij}(x_i - \mu_j)(x_i - \mu_j)</math> в вычислении ковариационных матриц <math>\Sigma_{j}</math> в M-шаге вычисляются независимо (<math>i = \overline{1, n}, \, j = \overline{1, k}</math>). Таким образом, при условии наличия неограниченного числа процессоров эти операции можно распараллелить как по числу наблюдений, так и по числу кластеров. Однако в отсутствии неограниченного числа процессорв целесообразно проводить распараллеливание по числу наблюдений, потому что, как правило, число наблюдений значительно больше числа кластеров  (<math>n \gg q</math>).
 +
 +
При классификации по ширине ЯПФ его сложность будет <math>O(nk)</math>. Сложность ЯПФ при классификации по высоте зависит от способа вычисления обратных матриц и определителей ковариационных матриц, а также от способа суммирования при вычислении скрытых переменных и обновлении вектора параметров. Предположим, что размерность элементов выборки <math>q</math> и число кластеров <math>k</math> значительно меньше числа элементов выборки <math>n</math> (<math>k \ll n, q \ll n</math>). Тогда вычисление обратных матриц и определителей ковариационных матриц можно не учитывать. Пусть для вычисления сумм используется пирамида сдваиванием. Тогда сложность суммирования при вычислении скрытых переменных равна <math>O(\log{k})</math>, а при обновлении вектора параметров — <math>O(\log{n})</math>. Следовательно, при классификации по высоте ЯПФ его сложность будет <math>O(\log{nk})</math>.
  
 
=== Входные и выходные данные алгоритма ===
 
=== Входные и выходные данные алгоритма ===
Строка 220: Строка 232:
 
=== Свойства алгоритма ===
 
=== Свойства алгоритма ===
 
Алгоритм является ''недетерминированным''. Число итераций до сходимости метода зависит от начального приближения и выбранного критерия останова.
 
Алгоритм является ''недетерминированным''. Число итераций до сходимости метода зависит от начального приближения и выбранного критерия останова.
 +
 +
Предположим, что размерность элементов выборки <math>q</math> и число кластеров <math>k</math> значительно меньше числа элементов выборки <math>n</math> (<math>k \ll n, q \ll n</math>), а число итераций до сходимости метода равно <math>m</math>. Тогда вычислительная мощность алгоритма, как отношение числа операций к суммарному объему входных и выходных данных, имеет порядок <math>O\bigg(\frac{kq^2}{q + k}m\bigg)</math>.
 
==== Достоинства и недостатки алгоритма ====
 
==== Достоинства и недостатки алгоритма ====
 
Среди достоинств EM-алгоритма можно выделить следующие:
 
Среди достоинств EM-алгоритма можно выделить следующие:
Строка 234: Строка 248:
 
== Программная реализация алгоритма ==
 
== Программная реализация алгоритма ==
  
=== Масштабируемость реализации алгоритма ===
+
=== Особенности реализации последовательного алгоритма ===
 +
 
 +
=== Локальность данных и вычислений ===
 +
 
 +
=== Возможные способы и особенности параллельной реализации алгоритма ===
 +
 
 +
Для исследования масштабируемости EM-алгоритм кластеризации был реализован с помощью технологии MPI. Была реализована версия алгоритма с диагональными матрицами ковариации. В качестве критерия останова было выбрано достижение максимального числа итераций (соответственно, вычисление логарифмического правдоподобия на каждой итерации не проводится).
 +
 
 +
Алгоритм параллельной реализации:
 +
# Распределение элементов выборки. Процесс 0 распределяет поровну элементы выборки по остальным процессам.
 +
# Инициализация параметров смеси гауссовых распределений. Процесс 0 инициализирует параметры и отправляет их остальным процессам. <math>\mu_j</math> инициализируются случайными значениями в отрезке  <math>[-1, 1]</math>, ковариационные матрицы <math>\Sigma_j</math> определяются как единичные, веса <math>w_i</math> кластеров задаются равными <math>\frac{1}{k}</math>.
 +
# Выполнение E-шага и M-шага алгоритма до достижения максимального числа итераций.
 +
## E-шаг. Каждый процесс вычисляет значения скрытых переменных <math>y_{ij}</math> по своему набору элементов выборки.
 +
## M-шаг.
 +
### Каждый процесс вычисляет суммы значений скрытых переменных <math>\sum_{i = 1}^{n\_proc_i}y_{ij}</math> для каждого кластера <math>j = \overline{1, k}</math> по своему набору элементов выборки и отправляет полученные суммы процессу 0. Процесс 0 вычисляет значения <math>n_j</math> по полученным суммам и отправляет их остальным процессам. Каждый процесс обновляет значения <math>w_j</math> по полученным значениям <math>n_j</math>.
 +
### Каждый процесс вычисляет суммы <math>\sum_{i = 1}^{n\_proc_i}y_{ij}x_i</math> для каждого кластера <math>j = \overline{1, k}</math> по своему набору элементов выборки и отправляет полученные суммы процессу 0. Процесс 0 вычисляет суммы <math>\sum_{i = 1}^{n}y_{ij}x_i</math> для каждого кластера <math>j = \overline{1, k}</math> по всем элементам выборки и отправляет их остальным процессам. С помощью полученных сумм каждый процесс вычисляет значения математических ожиданий <math>\mu_j</math>.
 +
### Каждый процесс вычисляет суммы <math>sum_l = \sum_{i = 1}^{n\_proc_i}y_{ij}(x_{il} - \mu_{jl}^{new})^2 \, (l = \overline{1, q})</math> для каждого кластера <math>j = \overline{1, k}</math> по своему набору элементов выборки и отправляет полученные суммы процессу 0. Процесс 0 вычисляет суммы <math>sum_l = \sum_{i = 1}^{n}y_{ij}(x_{il} - \mu_{jl}^{new})^2 \, (l = \overline{1, q})</math> для каждого кластера <math>j = \overline{1, k}</math> по всем элементам выборки и отправляет их остальным процессам. С помощью полученных сумм каждый процесс вычисляет значения диагональных элементов матриц ковариации <math>\Sigma_j</math>.
 +
 
 +
=== Масштабируемость алгоритма и его реализации ===
 +
 
 +
==== Масштабируемость алгоритма ====
 +
 
 +
==== Масштабируемость реализации алгоритма ====
 +
 
 +
Распараллеливание было проведено по числу элементов выборки, поэтому и исследование масштабируемости было проведено для числа элементов выборки. Исследование проводилось на суперкомпьютере "Ломоносов"<ref name="Lom">Воеводин Вл., Жуматий С., Соболев С., Антонов А., Брызгалов П., Никитенко Д., Стефанов К., Воеводин Вад. Практика суперкомпьютера «Ломоносов» // Открытые системы, 2012, N 7, С. 36-39.</ref> [http://parallel.ru/cluster Суперкомпьютерного комплекса Московского университета].
 +
 
 +
Набор и границы значений изменяемых параметров запуска реализации алгоритма:
 +
* число процессов [8 : 128] с шагом 8;
 +
* число элементов выборки [1000000 : 20000000] с шагом 1000000;
 +
* размерность элементов выборки 4;
 +
* число кластеров 4;
 +
* число итераций 50.
 +
 
 +
Данные, необходимые для тестирования алгоритма, генерируются в самой программе. Для того, чтобы данные были одинаковыми при каждом запуске, ядро рандомизации было зафиксировано.
 +
 
 +
Для каждого набора параметров было проведено по 3 запуска для избежания выбросов. Для каждого запуска измерялось среднее время выполнения одной итерации алгоритма. На следующем рисунке приведен график производительности данной реализации EM-алгоритма кластеризации в зависимости от изменяемых параметров запуска.
 +
 
 +
[[Файл:Em-performance-avg-Lavrushkin.png|thumb|1650px|center|Рисунок 8. Параллельная реализация EM-алгоритма кластеризации. Изменение производительности в зависимости от числа процессов и числа элементов выборки.]]
 +
 
 +
Данный алгоритм хорошо масштабируется по числу элементов выборки. По графику видно, что в целом производительность увеличивается при увеличении числа процессов.
 +
 
 +
[https://github.com/HighVoltageRocknRoll/EM-clusterization_scalability Исследованная параллельная MPI-реализация на языке C]
 +
 
 +
=== Динамические характеристики и эффективность реализации алгоритма ===
 +
 
 +
=== Выводы для классов архитектур ===
  
 
=== Существующие реализации алгоритма ===
 
=== Существующие реализации алгоритма ===
# [http://docs.opencv.org/2.4/modules/ml/doc/expectation_maximization.html OpenCV]
+
# [http://docs.opencv.org/2.4/modules/ml/doc/expectation_maximization.html OpenCV];
# [http://arma.sourceforge.net/docs.html#gmm_diag Armadillo] - OpenMP реализация с предположением диагональности ковариационных матриц
+
# [http://arma.sourceforge.net/docs.html#gmm_diag Armadillo] OpenMP реализация с предположением диагональности ковариационных матриц;
# [http://www.mlpack.org/docs/mlpack-git/doxygen.php?doc=classmlpack_1_1gmm_1_1GMM.html MLPack]
+
# [http://www.mlpack.org/docs/mlpack-git/doxygen.php?doc=classmlpack_1_1gmm_1_1GMM.html MLPack];
# [http://scikit-learn.org/stable/modules/mixture.html scikit-learn]
+
# [http://scikit-learn.org/stable/modules/mixture.html scikit-learn];
# [http://www.mathworks.com/help/stats/clustering-using-gaussian-mixture-models.html Matlab]
+
# [http://www.mathworks.com/help/stats/clustering-using-gaussian-mixture-models.html Matlab];
# [https://github.com/Corv/CUDA-GMM-MultiGPU CUDA реализация №1]
+
# [https://github.com/Corv/CUDA-GMM-MultiGPU CUDA реализация №1];
# [http://andrewharp.com/gmmcuda CUDA реализация №2]
+
# [http://andrewharp.com/gmmcuda CUDA реализация №2].
  
 
== Литература ==
 
== Литература ==
 
<references \>
 
<references \>

Текущая версия на 15:49, 28 ноября 2016

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

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



EM-алгоритм кластеризации (1 итерация)
Последовательный алгоритм
Последовательная сложность [math]O(nkq^2)[/math]
Объём входных данных [math]nq + 1[/math]
Объём выходных данных [math]k(n + q^2 + q + 1)[/math]
Параллельный алгоритм
Высота ярусно-параллельной формы [math]O(\log{nk})[/math]
Ширина ярусно-параллельной формы [math]O(nk)[/math]

Автор статьи: Сергей Лаврушкин (группа 620)

Содержание

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

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

1.1.1 EM-алгоритм

EM-алгоритм (англ. expectation-maximization) — алгоритм, используемый в математической статистике для нахождения оценок максимального правдоподобия параметров вероятностных моделей, в случае, когда модель зависит от некоторых скрытых переменных. Каждая итерация алгоритма состоит из двух шагов. На E-шаге (expectation) вычисляется ожидаемое значение функции правдоподобия, при этом скрытые переменные рассматриваются как наблюдаемые. На M-шаге (maximization) вычисляется оценка максимального правдоподобия, таким образом увеличивается ожидаемое правдоподобие, вычисляемое на E-шаге. Затем это значение используется для E-шага на следующей итерации. Алгоритм выполняется до сходимости.[1]

Как правило, ЕМ-алгоритм применяется для решения задач двух типов:

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

Алгоритм основан на методике итеративного вычисления оценок максимального правдоподобия, предложенной в 1977 г. (A. P. Demster, N. M. Laird, D. B. Rubin. Maximum Likelihood from Incomplete Data via the EM Algorithm).[2]

1.1.2 Кластеризация с помощью EM-алгоритма

Кластеризация — задача разбиения заданной выборки объектов (ситуаций) на непересекающиеся подмножества, называемые кластерами, так, чтобы каждый кластер состоял из схожих объектов, а объекты разных кластеров существенно отличались [3]. В основе идеи применения EM-алгоритма для кластеризации лежит предположение, что исследуемое множество данных может быть смоделировано с помощью линейной комбинации многомерных нормальных распределений, а целью является оценка параметров распределения, которые максимизируют логарифмическую функцию правдоподобия, используемую в качестве меры качества модели. К оцениваемым параметрам относятся математические ожидания и ковариационные матрицы. Предполагается, что любое наблюдение принадлежит ко всем кластерам, но с разной вероятностью. Задача заключается в "подгонке" распределений смеси к данным, а затем в определении вероятностей принадлежности наблюдения к каждому кластеру. Наблюдение должно быть отнесено к тому кластеру, для которого данная вероятность выше.

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

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

Исходные данные:

  • [math]k[/math] — число кластеров,
  • [math]X =\{x_{1}, x_{2}, ..., x_{n}\}[/math] — множество из [math]n[/math] наблюдений размерности [math]q[/math]: [math]x_{i} = (x_{i1}, ..., x_{iq})[/math].

Вычисляемые данные:

  • [math](w, \theta) = (w_1, ..., w_k; \mu_1, ..., \mu_k; \Sigma_1, ..., \Sigma_k)[/math] — параметры смеси гауссовых распределений,
  • [math]Y[/math] — матрица размера [math]n \times k[/math] с вероятностями членства в кластерах.

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

Пусть [math]w_1, ..., w_k[/math] — априорные вероятности кластеров, [math]p_1(x), ..., p_k(x)[/math] — плотности распределения кластеров. Тогда плотность распределения вектора признаков [math]x[/math] сразу по всем кластерам равна:

[math] \begin{align} p(x) = \sum_{j=1}^{k}w_jp_j(x) \end{align} [/math]

Необходимо на основе выборки оценить параметры модели [math]w_1, ..., w_k, p_1(x), ..., p_k(x)[/math]. Это позволит оценивать вероятность принадлежности к кластеру и, таким образом, решить задачу кластеризации. Такая задача называется задачей разделения смеси распределений:

[math] \begin{align} p(x) = \sum_{j=1}^{k}w_jp_j(x),\quad p_j(x) = \phi(\theta_j; x), \end{align} [/math]

где [math]\theta_j[/math] — параметры распределения [math]p_j(x)[/math]. В случае смеси гауссовых распределений в [math]q[/math]-мерном пространстве признаков [math]\phi(\theta_j; x)[/math] имеют следующий вид:

[math] \begin{align} \phi(\theta_j; x) = \frac{1}{\Big(2\pi\Big)^{\frac{q}{2}}\sqrt{|\Sigma_j|}}\exp{\bigg\{-\frac{\delta^2}{2}\bigg\} }, \end{align} [/math]

где:

  • [math]\Sigma_j[/math] — ковариационная матрица размером [math]q \times q[/math],
  • [math]\mu_j[/math][math]q[/math]-мерный вектор математических ожиданий,
  • [math]\delta^2 = \Big( x - \mu_j \Big)^T \Sigma_j^{-1} \Big( x - \mu_j \Big)[/math] — квадратичное расстояние Махаланобиса.

В случае смеси гауссовых распределений получаем, что [math]\theta_j = (\mu_j, \Sigma_j)[/math]. То есть для решения задачи разделения смеси гауссовых распределений необходимо оценить вектор параметров [math](w, \theta) = (w_1, ..., w_k; \mu_1, ..., \mu_k; \Sigma_1, ..., \Sigma_k)[/math]. Согласно принципу максимизации правдоподобия:

[math] \begin{align} w, \theta = \underset{w, \theta}{\operatorname{argmax}} \sum_{i=1}^{n}\ln{p(x_i)} = \underset{w, \theta}{\operatorname{argmax}} \sum_{i=1}^{n}\ln{\sum_{j=1}^{k}w_jp_j(x_i)} \end{align} [/math]

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

1.2.2 EM-алгоритм

Оптимальные параметры задачи разделения смеси гауссовых распределений отыскиваются последовательно с помощью итерационного EM-алгоритма. Основная идея – вводится вспомогательный вектор скрытых переменных. Это позволяет свести сложную оптимизационную задачу к последовательности итераций по пересчету коэффициентов (скрытых переменных по текущему приближению вектора параметров - E-шаг) и максимизации правдоподобия (с целью найти следующее приближение вектора - М-шаг).[6]

EM-алгоритм заключается в следующем:

  1. В начале работы алгоритма задаются параметры начального приближения. Наиболее общим способом инициализации является присвоение элементам матрицы математических ожиданий случайных значений [math] \mu_j \leftarrow Random [/math], начальные ковариационные матрицы определяются как единичные [math]\Sigma_j \leftarrow I[/math], веса кластеров задаются одинаковыми [math]w_i \leftarrow \frac{1}{k}[/math]. Также в качестве начальных параметров можно использовать результат работы алгоритма K-means. (Данная эвристика применяется, так как K-means требуется намного меньше итераций до достижения стабилизации, в то время как каждый шаг EM требует больших вычислительных затрат).
  2. Далее итеративно выполняется следующая пара процедур:
    • E-шаг: используя текущее значение вектора параметров [math](w, \theta) = (w_1, ..., w_k; \mu_1, ..., \mu_k; \Sigma_1, ..., \Sigma_k)[/math], вычисляются значения скрытых переменных [math]y_{ij}[/math]:
      [math] \begin{align} & y_{ij} = \frac{w_j p_j(x_i)}{\sum_{l=1}^{k}w_lp_l(x_i)} \end{align} [/math]
    • М-шаг: переоценка вектора параметров, используя текущие значения скрытых переменных:
      [math] \begin{align} & n_j = \sum_{i = 1}^{n}y_{ij}, \\ & \mu_j^{new} = \frac{1}{n_j} \sum_{i = 1}^{n}y_{ij}x_i, \\ & \Sigma_j^{new} = \frac{1}{n_j} \sum_{i = 1}^{n}y_{ij}(x_i - \mu_j^{new})(x_i - \mu_j^{new})^T,\\ & w_j^{new} = \frac{n_j}{n}. \end{align} [/math]

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

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

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

  • Вычисление скрытых переменных [math]y_{ij}[/math] на E-шаге (всего их [math]nk[/math]).
  • Переоценка вектора параметров [math](w_1, ..., w_k; \mu_1, ..., \mu_k; \Sigma_1, ..., \Sigma_k)[/math] на M-шаге, используя текущие значения скрытых переменных [math]y_{ij}[/math].

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

Как записано в описании ядра алгоритма, основную часть метода составляют множественные вычисления скрытых переменных [math]y_{ij}[/math] на E-шаге и обновление вектора параметров [math](w_1, ..., w_k; \mu_1, ..., \mu_k; \Sigma_1, ..., \Sigma_k)[/math] на M-шаге. Также стоит выделить вычисление определителя и обратной матрицы ковариационных матриц [math]\Sigma_j[/math] в начале E-шага.

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

Последовательная реализация алгоритма EM может быть проиллюстрирована с помощью следующего псевдокода:

  1. Инициализация: установка начальных значений [math]\mu, \, \Sigma_j \, (j = \overline{1, k}), \, w[/math].
  2. Пока изменение логарифмического правдоподобия [math]\Delta llh \geq \epsilon[/math] и не достигнуто максимальное число итераций [math]m[/math], выполнять шаги E и M.

Шаг E

def E-step:
    Инициализировать нулевыми значениями [math]\mu^{'}, \, \Sigma_j^{'} \, (j = \overline{1, k}), \, w^{'}, \, llh[/math]
    for [math]j = \overline{1, k}[/math]:
        Вычислить определитель и обратную матрицу матрицы [math]\Sigma_j[/math]
    for [math]i = \overline{1, n}[/math]:
        [math]sump_i = 0[/math]
        for [math]j = \overline{1, k}[/math]:
            [math]\delta_{ij}^2 = \big( x_i - \mu_j \big)^T \Sigma_j^{-1} \big( x_i - \mu_j \big)[/math]
            [math]g_{ij} = \frac{w_j}{{(2 \pi)}^\frac{q}{2}\sqrt{|\Sigma_j|}}\exp{\big\{-\frac{\delta_{ij}^2}{2}\big\} }[/math]
            [math]sump_i \mathrel{+}= g_{ij}[/math]
        [math]y_i = \frac{g_i}{sump_i}[/math]
        [math]llh \mathrel{+}= \ln(sump_i)[/math]
        [math]\mu^{'} \mathrel{+}= x_iy_i^T[/math]
        [math]w^{'} \mathrel{+}= y_i[/math]

Шаг M

def M-step:
    for [math]j = \overline{1, k}[/math]:
        [math]\mu_j = \frac{\mu_j^{'}}{w_j^{'}}[/math]
        for [math]i = \overline{1, n}[/math]:
            [math]\Sigma_j^{'} \mathrel{+}= (x_i - \mu_j)y_{ij}(x_i - \mu_j)^T[/math]
        [math]\Sigma_j = \frac{\Sigma_j^{'}}{w_j^{'}}[/math]
    [math]w = \frac{w^{'}}{n}[/math]

Логаририфмическое правдоподобие вычисляется как:

[math] \begin{align} llh=\sum_{i=1}^{n}\ln{(sump_i)} \end{align} [/math]

Матрицы [math]\mu^{'}, \, \Sigma_j^{'} \, (j = \overline{1, k}), \, w^{'}[/math] являются временными и используются только для вычислений.

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

Последовательная сложность одной итерации EM-алгоритма:

  1. E-шаг:
    В случае недиагональных ковариационных матриц сложность вычислений:
    • Определителей ковариационных матриц [math]|\Sigma_j|[/math] ([math]j = \overline{1, k}[/math]) методом Гаусса: [math]O(kq^3)[/math],
    • Матриц, обратных ковариационным матрицам [math]\Sigma_j^{-1}[/math] ([math]j = \overline{1, k}[/math]) методом Гаусса-Жордана: [math]O(kq^3)[/math],
    • Скрытых переменных [math]y_{ij}[/math] ([math]i = \overline{1, n}, j = \overline{1, k}[/math]): [math]O(nkq^2)[/math].
    В случае диагональных ковариационных матриц сложность вычислений:
    • Определителей ковариационных матриц [math]|\Sigma_j|[/math] ([math]j = \overline{1, k}[/math]): [math]O(kq)[/math],
    • Матриц, обратных ковариационным матрицам [math]\Sigma_j^{-1}[/math] ([math]j = \overline{1, k}[/math]): [math]O(kq)[/math],
    • Скрытых переменных [math]y_{ij}[/math] ([math]i = \overline{1, n}, j = \overline{1, k}[/math]): [math]O(nkq)[/math].
    Получаем, что последовательная сложность E-шага в случае недиагональных ковариационных матриц имеет порядок [math]O(kq^3 + nkq^2)[/math], а в случае диагональных ковариационных матриц — [math]O(nkq)[/math].
  2. М-шаг:
    В случае недиагональных ковариационных матриц сложность вычислений:
    • Обновление математических ожиданий [math]\mu_j[/math] ([math]j = \overline{1, k}[/math]): [math]O(nkq)[/math],
    • Обновление ковариационных матриц [math]\Sigma_j[/math] ([math]j = \overline{1, k}[/math]): [math]O(nkq^2)[/math],
    • Обновление весов кластеров [math]w_j[/math] ([math]j = \overline{1, k}[/math]): [math]O(nk)[/math].
    В случае диагональных ковариационных матриц сложность вычислений:
    • Обновление математических ожиданий [math]\mu_j[/math] ([math]j = \overline{1, k}[/math]): [math]O(nkq)[/math],
    • Обновление ковариационных матриц [math]\Sigma_j[/math] ([math]j = \overline{1, k}[/math]): [math]O(nkq)[/math],
    • Обновление весов кластеров [math]w_j[/math] ([math]j = \overline{1, k}[/math]): [math]O(nk)[/math].
    Получаем, что последовательная сложность М-шага в случае недиагональных ковариационных матриц имеет порядок [math]O(nkq^2)[/math], а в случае диагональных ковариационных матриц — [math]O(nkq)[/math].

Таким образом, последовательная сложность итерации EM-алгоритма в случае недиагональных ковариационных матриц имеет порядок [math]O(kq^3 + nkq^2)[/math], а в случае диагональных ковариационных матриц — [math]O(nkq)[/math]. При этом, обычно размерность элементов выборки [math]q[/math], также как и число кластеров [math]k[/math], значительно меньше числа элементов выборки [math]n[/math] ([math]k \ll n, q \ll n[/math]), поэтому в случае недиагональных ковариационных матриц сложность может быть записана в виде [math]O(nkq^2)[/math].

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

Граф алгоритма в терминах E-шагов и M-шагов представлен на рис. 1. Каждая итерация алгоритма зависит от результата работы предыдущей итерации, M-шаг зависит от результата E-шага.

Рисунок 1. Граф EM-алгоритма. [math]k[/math] — число кластеров; [math]X[/math] — множество наблюдений; [math]Y[/math] — матрица с вероятностями членства в кластерах.

Рассмотрим по отдельности информационные графы E-шага и M-шага.

1.7.1 E-шаг

Опишем граф E-шага как аналитически, так и в виде рисунка.

Граф алгоритма состоит из четырех групп вершин:

  • Первая группа вершин соответствует вычислению обратных матриц и определителей ковариационных матриц. Число вершин: [math]k[/math]. Аргументом этой операции является ковариационная матрица [math]\Sigma_j[/math]. Результат срабатывания операции является промежуточным данным алгоритма.
  • Вторая группа вершин соответствует вычислению [math]g_{ij} = \frac{w_j}{{2 \pi}^\frac{q}{2}|\Sigma_j|}\exp{\big\{-\frac{1}{2} \big( x_i - \mu_j \big)^T \Sigma_j^{-1} \big( x_i - \mu_j \big) \big\}}[/math]. Число вершин: [math]nk[/math]. Аргументами этой операции являются вес кластера [math]w_j[/math], математическое ожидание [math]\mu_j[/math], ковариационная матрица [math]\Sigma_j[/math], элемент выборки [math]x_i[/math]. Результат срабатывания операции является промежуточным данным алгоритма.
  • Третья группа вершин соответствует вычислению нормировочных коэффициентов [math]\sum_{l=1}^{k}g_{ij}[/math]. Число вершин: [math]n[/math]. Аргументами этой операции являются значения [math]g_{ij}, \, j = \overline{1, k}[/math]. Результат срабатывания операции является промежуточным данным алгоритма.
  • Четвертая группа вершин соответствует вычислению скрытых переменных [math]y_{ij}[/math]. Число вершин: [math]nk[/math]. Аргументами этой операции являются значение [math]g_{ij}[/math] и нормировочный коэффициент. Результат срабатывания операции является выходным данным алгоритма.

Описанный граф представлен на рис.2 и рис.3. Здесь вершины одной группы расположены в одной строке. Вершины-операции обозначены синим цветом, вершины-данные — оранжевым.

Рисунок 2. Граф E-шага.
Рисунок 3. Граф E-шага с с отображением входных и выходных данных.

1.7.2 M-шаг

Опишем граф M-шага как аналитически, так и в виде рисунка.

Граф алгоритма состоит из четырех групп вершин:

  • Первая группа вершин соответствует вычислению нормировочных коэффициентов [math]n_j[/math]. Число вершин: [math]k[/math]. Аргументами этой операции являются скрытые переменные [math]y_{ij}, \, i = \overline{1, n}[/math]. Результат срабатывания операции является промежуточным данным алгоритма.
  • Вторая группа вершин соответствует вычислению математических ожиданий [math]\mu_j[/math]. Число вершин: [math]k[/math]. Аргументами этой операции являются множество наблюдений [math]X[/math], скрытые переменные [math]y_{ij}, \, i = \overline{1, n}[/math], нормировочный коэффициент [math]n_j[/math]. Результат срабатывания операции является выходным данным алгоритма.
  • Третья группа вершин соответствует вычислению ковариационных матриц [math]\Sigma_j[/math]. Число вершин: [math]k[/math]. Аргументами этой операции являются множество наблюдений [math]X[/math], скрытые переменные [math]y_{ij}, \, i = \overline{1, n}[/math], нормировочный коэффициент [math]n_j[/math], математическое ожидание [math]\mu_j[/math]. Результат срабатывания операции является выходным данным алгоритма.
  • Четвертая группа вершин соответствует вычислению весов кластеров [math]w_{j}[/math]. Число вершин: [math]k[/math]. Аргументом этой операции является нормировочный коэффициент [math]n_j[/math]. Результат срабатывания операции является выходным данным алгоритма.

Описанный граф представлен на рис.4 и рис.5. Здесь вершины одной группы расположены в одной строке. Вершины-операции обозначены синим цветом, вершины-данные — оранжевым.

Рисунок 4. Граф M-шага.
Рисунок 5. Граф M-шага с отображением входных и выходных данных.

Графы вычисления математических ожиданий [math]\mu_j[/math] и ковариационных матриц [math]\Sigma_j[/math] представлены на рис. 6 и рис. 7.

Рисунок 6. Граф вычисления математического ожидания [math]\mu_j[/math].
Рисунок 7. Граф вычисления ковариационной матрицы [math]\Sigma_j[/math].

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

По графам из пункта 1.7 видно, что скрытые переменные [math]y_{ij}[/math] в E-шаге, [math]y_{ij}x_i[/math] в вычислении математических ожиданий [math]\mu_{j}[/math] и [math]y_{ij}(x_i - \mu_j)(x_i - \mu_j)[/math] в вычислении ковариационных матриц [math]\Sigma_{j}[/math] в M-шаге вычисляются независимо ([math]i = \overline{1, n}, \, j = \overline{1, k}[/math]). Таким образом, при условии наличия неограниченного числа процессоров эти операции можно распараллелить как по числу наблюдений, так и по числу кластеров. Однако в отсутствии неограниченного числа процессорв целесообразно проводить распараллеливание по числу наблюдений, потому что, как правило, число наблюдений значительно больше числа кластеров ([math]n \gg q[/math]).

При классификации по ширине ЯПФ его сложность будет [math]O(nk)[/math]. Сложность ЯПФ при классификации по высоте зависит от способа вычисления обратных матриц и определителей ковариационных матриц, а также от способа суммирования при вычислении скрытых переменных и обновлении вектора параметров. Предположим, что размерность элементов выборки [math]q[/math] и число кластеров [math]k[/math] значительно меньше числа элементов выборки [math]n[/math] ([math]k \ll n, q \ll n[/math]). Тогда вычисление обратных матриц и определителей ковариационных матриц можно не учитывать. Пусть для вычисления сумм используется пирамида сдваиванием. Тогда сложность суммирования при вычислении скрытых переменных равна [math]O(\log{k})[/math], а при обновлении вектора параметров — [math]O(\log{n})[/math]. Следовательно, при классификации по высоте ЯПФ его сложность будет [math]O(\log{nk})[/math].

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

Входные данные: [math]k[/math] — число кластеров; [math]X =\{x_{1}, x_{2}, ..., x_{n}\}[/math] — множество из [math]n[/math] наблюдений [math]q[/math]-мерного пространства.

Объем входных данных: [math]nq + 1[/math].

Выходные данные: [math](w, \theta) = (w_1, ..., w_k; \mu_1, ..., \mu_k; \Sigma_1, ..., \Sigma_k)[/math] — параметры смеси гауссовых распределений; [math]Y[/math] — матрица размера [math]n \times k[/math] с вероятностями членства в кластерах.

Объем выходных данных: [math]k(n + q^2 + q + 1)[/math].

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

Алгоритм является недетерминированным. Число итераций до сходимости метода зависит от начального приближения и выбранного критерия останова.

Предположим, что размерность элементов выборки [math]q[/math] и число кластеров [math]k[/math] значительно меньше числа элементов выборки [math]n[/math] ([math]k \ll n, q \ll n[/math]), а число итераций до сходимости метода равно [math]m[/math]. Тогда вычислительная мощность алгоритма, как отношение числа операций к суммарному объему входных и выходных данных, имеет порядок [math]O\bigg(\frac{kq^2}{q + k}m\bigg)[/math].

1.10.1 Достоинства и недостатки алгоритма

Среди достоинств EM-алгоритма можно выделить следующие:

  • Мощная статистическая основа.
  • Линейное увеличение сложности при росте объема данных.
  • Устойчивость к шумам и пропускам в данных.
  • Возможность построения желаемого числа кластеров.
  • Быстрая сходимость при удачной инициализации.

К недостаткам EM-алгоритма можно отнести следующее:

  • Алгоритм неустойчив по начальным данным (то есть тем, которые инициализируют вектор параметров на первой итерации), так как он находит локальный экстремум, значение которого может оказаться гораздо ниже, чем глобальный максимум. В зависимости от выбора начального приближения алгоритм может сходиться к разным точкам. Также может сильно варьироваться скорость сходимости.
  • Алгоритм не позволяет определять количество [math]k[/math] компонент смеси. Эта величина является структурным параметром алгоритма.
  • Предположение о нормальности всех измерений данных не всегда выполняется.

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

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

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

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

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

Алгоритм параллельной реализации:

  1. Распределение элементов выборки. Процесс 0 распределяет поровну элементы выборки по остальным процессам.
  2. Инициализация параметров смеси гауссовых распределений. Процесс 0 инициализирует параметры и отправляет их остальным процессам. [math]\mu_j[/math] инициализируются случайными значениями в отрезке [math][-1, 1][/math], ковариационные матрицы [math]\Sigma_j[/math] определяются как единичные, веса [math]w_i[/math] кластеров задаются равными [math]\frac{1}{k}[/math].
  3. Выполнение E-шага и M-шага алгоритма до достижения максимального числа итераций.
    1. E-шаг. Каждый процесс вычисляет значения скрытых переменных [math]y_{ij}[/math] по своему набору элементов выборки.
    2. M-шаг.
      1. Каждый процесс вычисляет суммы значений скрытых переменных [math]\sum_{i = 1}^{n\_proc_i}y_{ij}[/math] для каждого кластера [math]j = \overline{1, k}[/math] по своему набору элементов выборки и отправляет полученные суммы процессу 0. Процесс 0 вычисляет значения [math]n_j[/math] по полученным суммам и отправляет их остальным процессам. Каждый процесс обновляет значения [math]w_j[/math] по полученным значениям [math]n_j[/math].
      2. Каждый процесс вычисляет суммы [math]\sum_{i = 1}^{n\_proc_i}y_{ij}x_i[/math] для каждого кластера [math]j = \overline{1, k}[/math] по своему набору элементов выборки и отправляет полученные суммы процессу 0. Процесс 0 вычисляет суммы [math]\sum_{i = 1}^{n}y_{ij}x_i[/math] для каждого кластера [math]j = \overline{1, k}[/math] по всем элементам выборки и отправляет их остальным процессам. С помощью полученных сумм каждый процесс вычисляет значения математических ожиданий [math]\mu_j[/math].
      3. Каждый процесс вычисляет суммы [math]sum_l = \sum_{i = 1}^{n\_proc_i}y_{ij}(x_{il} - \mu_{jl}^{new})^2 \, (l = \overline{1, q})[/math] для каждого кластера [math]j = \overline{1, k}[/math] по своему набору элементов выборки и отправляет полученные суммы процессу 0. Процесс 0 вычисляет суммы [math]sum_l = \sum_{i = 1}^{n}y_{ij}(x_{il} - \mu_{jl}^{new})^2 \, (l = \overline{1, q})[/math] для каждого кластера [math]j = \overline{1, k}[/math] по всем элементам выборки и отправляет их остальным процессам. С помощью полученных сумм каждый процесс вычисляет значения диагональных элементов матриц ковариации [math]\Sigma_j[/math].

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

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

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

Распараллеливание было проведено по числу элементов выборки, поэтому и исследование масштабируемости было проведено для числа элементов выборки. Исследование проводилось на суперкомпьютере "Ломоносов"[7] Суперкомпьютерного комплекса Московского университета.

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

  • число процессов [8 : 128] с шагом 8;
  • число элементов выборки [1000000 : 20000000] с шагом 1000000;
  • размерность элементов выборки 4;
  • число кластеров 4;
  • число итераций 50.

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

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

Рисунок 8. Параллельная реализация EM-алгоритма кластеризации. Изменение производительности в зависимости от числа процессов и числа элементов выборки.

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

Исследованная параллельная MPI-реализация на языке C

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

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

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

  1. OpenCV;
  2. Armadillo — OpenMP реализация с предположением диагональности ковариационных матриц;
  3. MLPack;
  4. scikit-learn;
  5. Matlab;
  6. CUDA реализация №1;
  7. CUDA реализация №2.

3 Литература

<references \>

  1. ЕМ-алгоритм, его модификации и обобщения http://www.machinelearning.ru/wiki/index.php?title=EM-алгоритм
  2. A. Dempster, N. Laird and D. Rubin. Maximum likelihood estimation from incomplete data. – Journal of the Royal Statistical Society, Series B, 1977, vol. 39, p. 1-38.
  3. Кластеризация http://www.machinelearning.ru/wiki/index.php?title=Кластеризация
  4. EМ — масштабируемый алгоритм кластеризации https://basegroup.ru/community/articles/em
  5. Поиск структуры в данных — Lecture 9 — Expectation Maximization (EM-алгоритм) https://www.coursera.org/learn/unsupervised-learning/lecture/2jhbI/expectation-maximization-em-alghoritm
  6. EM алгоритм (пример) http://www.machinelearning.ru/wiki/index.php?title=EM_алгоритм_(пример)
  7. Воеводин Вл., Жуматий С., Соболев С., Антонов А., Брызгалов П., Никитенко Д., Стефанов К., Воеводин Вад. Практика суперкомпьютера «Ломоносов» // Открытые системы, 2012, N 7, С. 36-39.