Участник:Noite/EM-алгоритм кластеризации: различия между версиями
Noite (обсуждение | вклад) |
Algoman (обсуждение | вклад) |
||
(не показано 148 промежуточных версий 4 участников) | |||
Строка 1: | Строка 1: | ||
{{algorithm | {{algorithm | ||
− | | name = EM-алгоритм | + | | name = EM-алгоритм кластеризации |
− | | serial_complexity = <math></math> | + | | serial_complexity = <math>O(iter * N * M^2 * K)</math> |
− | | | + | | input_data = <math>N * M + 1</math> |
− | | | + | | output_data = <math>N</math> |
− | | | + | | pf_height = <math>O(iter * M^2 * K)</math> |
− | | | + | | pf_width = <math>O(iter * N)</math> |
}} | }} | ||
Автор описания: [[Участник:Noite|Зинченко Д.А.]] | Автор описания: [[Участник:Noite|Зинченко Д.А.]] | ||
− | + | = Свойства и структура алгоритма = | |
− | + | == Общее описание алгоритма == | |
− | + | В задачах оптимизации EM-алгоритмом называют итеративную процедуру поиска численного значения экстремума какой-либо функции. В частности, в статистике этот алгоритм используется для оценки максимального правдоподобия. Впервые название EM-алгоритм было предложено в 1977 году<ref>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.</ref>, однако его идеи были описаны и раньше<ref name=korolev>В.Ю.Королёв. ЕМ-алгоритм, его модификации и их применение к задаче разделения смесей вероятностных распределений.</ref>. | |
− | В задачах оптимизации EM-алгоритмом называют итеративную процедуру поиска численного значения экстремума какой-либо функции. В частности, в статистике этот алгоритм используется для оценки максимального правдоподобия. Впервые название EM-алгоритм было предложено в 1977 году | ||
− | EM-алгоритм применяется для решения задач двух типов | + | EM-алгоритм применяется для решения задач двух типов <ref name=korolev/>: |
− | + | # Анализ неполных данных (данных с пропусками). | |
− | + | # Оценка максимального правдоподобия в случае, если функцию правдоподобия трудно исследовать аналитически. В этом случае введение набора скрытых переменных может существенно упростить задачу. | |
− | |||
Задача кластеризации относится к задачам второго типа. В данной статье рассматривается следующая постановка задачи кластеризации: | Задача кластеризации относится к задачам второго типа. В данной статье рассматривается следующая постановка задачи кластеризации: | ||
− | <source>Дано множество объектов, каждый из которых представляет собой точку в | + | <source>Дано множество из N объектов, каждый из которых представляет собой точку в M-мерном метрическом пространстве. Каждому измерению пространства соответствует некоторое свойство объекта. Необходимо разбить это множество на K подмножеств (кластеров) так, чтобы, чтобы элементы каждого подмножества существенно отличались по свойствам от элементов других подмножеств.</source> |
Использование алгоритма EM подразумевает следующее ''предположение об объектах'': | Использование алгоритма EM подразумевает следующее ''предположение об объектах'': | ||
Строка 33: | Строка 31: | ||
Итерация алгоритма состоит из двух последовательных шагов. | Итерация алгоритма состоит из двух последовательных шагов. | ||
− | + | # (Expectation) Вычисляются новые ожидаемые значения скрытых переменных. | |
− | + | # (Maximization) Решается задача максимизации правдоподобия: По текущим значениям скрытых переменных обновляются параметры распределений для кластеров: математическое ожидание, дисперсия и вероятность появления объектов из кластера. | |
− | |||
− | + | == Математическое описание алгоритма == | |
''Используемые обозначения:'' | ''Используемые обозначения:'' | ||
Строка 45: | Строка 42: | ||
<math>\gamma_{n,k}</math> - значения скрытых переменных | <math>\gamma_{n,k}</math> - значения скрытых переменных | ||
− | <math>N_{k}</math> - количество элементов отнесенных к | + | <math>N_{k}</math> - количество элементов отнесенных к кластеру <math>k</math> |
− | <math>x_{1},..., | + | <math>x_{1},...,x_N</math> - объекты |
''Параметры распределений:'' | ''Параметры распределений:'' | ||
− | <math>w_{1},...,w_{ | + | <math>w_{1},...,w_{K}</math> - вероятности появления объектов из кластеров (веса кластеров), <math>\sum_{k=1}^Kw_{k}=1</math> |
− | <math>\mu_{1},...,\mu_{ | + | <math>\mu_{1},...,\mu_{K}</math> - центры гауссиан кластеров |
− | <math>\Sigma_{1},...,\Sigma_{ | + | <math>\Sigma_{1},...,\Sigma_{K}</math> - матрицы ковариации гауссиан кластеров, каждая матрица имеет размерность <math>M</math>x<math>M</math> |
''Смесь распределений:'' | ''Смесь распределений:'' | ||
Строка 70: | Строка 67: | ||
<center><math>\Sigma_{j,j}^{k}=\frac{1}{NK}\sum_{i=1}^N(x_{i,j} - \mu_{k,j})^2, k=1..K, j=1..M</math></center> | <center><math>\Sigma_{j,j}^{k}=\frac{1}{NK}\sum_{i=1}^N(x_{i,j} - \mu_{k,j})^2, k=1..K, j=1..M</math></center> | ||
− | ''Поскольку от начального приближения может сильно зависеть результат работы алгоритма, вместо инициализации кластеров случайными центрами, равными весами и диагональными матрицами ковариации, как в приведенных формулах, часто применяют другие способы (подробнее - в | + | ''Поскольку от начального приближения может сильно зависеть результат работы алгоритма, вместо инициализации кластеров случайными центрами, равными весами и диагональными матрицами ковариации, как в приведенных формулах, часто применяют другие способы (подробнее - в [[#Особенности реализации последовательного алгоритма | особенностях реализации]])'' |
* Далее итеративно выполняется следующая пара процедур: | * Далее итеративно выполняется следующая пара процедур: | ||
− | ** E-шаг: используя текущее значение параметров <math>w_{1},...,w_{k};\mu_{1},...,\mu_{k};\ | + | ** E-шаг: используя текущее значение параметров <math>w_{1},...,w_{k};\mu_{1},...,\mu_{k};\Sigma^{1},...,\Sigma^{k}</math>, вычисляем значение вектора скрытых переменных <math>\gamma</math>:<center><math>\gamma_{n,k}=\frac{w_{k}p(x_{n}|\mu_{k},\Sigma^{k})}{\sum_{j=1}^K w_{j}p(x_{n}|\mu_{j},\Sigma^{j})}, </math> знаменатель дроби - нормировочный коэффициент</center><center><math>p(x|\mu,\Sigma)=\frac{1}{{(2\pi)}^{\frac{N}{2}}\sqrt{|\Sigma|}}exp\biggl\{-\frac{1}{2}{(x-\mu)}^T{\Sigma}^{-1}(x-\mu)\biggr\}</math> - плотность <math>N</math>-мерного нормального распределения</center> |
− | ** М-шаг: переоценка параметров по текущим значениям скрытых переменных | + | ** М-шаг: переоценка параметров по текущим значениям скрытых переменных для <math>k = 1..K</math>: |
<center><math>\mu^{new}_{k}=\frac{1}{N_{k}}\sum_{n=1}^N \gamma_{n,k}x_{n}</math>,</center> | <center><math>\mu^{new}_{k}=\frac{1}{N_{k}}\sum_{n=1}^N \gamma_{n,k}x_{n}</math>,</center> | ||
Строка 85: | Строка 82: | ||
*Изменение логарифмического правдоподобия меньше заданной константы. | *Изменение логарифмического правдоподобия меньше заданной константы. | ||
− | + | == Вычислительное ядро алгоритма == | |
Вычислительное ядро алгоритма состоит из двух шагов, повторяющихся каждую итерацию: | Вычислительное ядро алгоритма состоит из двух шагов, повторяющихся каждую итерацию: | ||
*Обновление скрытых переменных на основе текущих параметров распределений | *Обновление скрытых переменных на основе текущих параметров распределений | ||
*Обновление параметров кластеров: весов, центров и матриц ковариации | *Обновление параметров кластеров: весов, центров и матриц ковариации | ||
− | + | == Макроструктура алгоритма == | |
− | Как записано в описании вычислительного ядра алгоритма, основную часть алгоритма составляют итеративно проводимые макрооперации: Expectation шаг - вычисление новых значений скрытых переменных, Maximization шаг - обновление параметров кластеров. | + | Как записано в [[#Вычислительное ядро алгоритма|описании вычислительного ядра алгоритма]], основную часть алгоритма составляют итеративно проводимые макрооперации: Expectation шаг - вычисление новых значений скрытых переменных, Maximization шаг - обновление параметров кластеров. Макроструктура алгоритма представлена на рисунке 1: |
+ | [[file:EM_macro.png|thumb|center|600px|Рисунок 1: Макроструктура EM-алгоритма]] | ||
+ | |||
+ | == Схема реализации последовательного алгоритма == | ||
+ | Примерная схема реализации EM алгоритма с предварительной обработкой матриц ковариации: | ||
+ | <source> | ||
+ | Инициализировать веса и центры кластеров; | ||
+ | Инициализировать матрицы ковариации; | ||
− | + | // Основной цикл | |
+ | Пока условие остановки не выполнено | ||
+ | // E-шаг | ||
+ | В цикле для каждого кластера | ||
+ | Вычислить определитель матрицы ковариации Sigma; | ||
+ | Вычислить обратную матрицу матрице Sigma; | ||
+ | В цикле для каждого объекта | ||
+ | В цикле для каждого кластера | ||
+ | Вычислить и запомнить p(x | mu, Sigma) для текущего объекта и кластера; | ||
+ | Просуммировать вычисленные p(x | mu, Sigma); | ||
+ | В цикле для каждого кластера | ||
+ | Обновить скрытые переменные gamma; | ||
+ | // M-шаг | ||
+ | Обновить веса w; | ||
+ | Обновить центры кластеров mu; | ||
+ | Обновить матрицы ковариации Sigma; | ||
+ | </source> | ||
− | + | == Последовательная сложность алгоритма == | |
Количество итераций алгоритма зависит от условия остановки, поэтому будем рассматривать сложность одной итерации. | Количество итераций алгоритма зависит от условия остановки, поэтому будем рассматривать сложность одной итерации. | ||
Строка 101: | Строка 121: | ||
'''Шаг E''' | '''Шаг E''' | ||
− | Сложность вычисления определителя произвольной матрицы ковариации <math>\Sigma</math> методом Гаусса - <math>O( | + | Сложность вычисления определителя произвольной матрицы ковариации <math>\Sigma</math> методом Гаусса - <math>O(M^3)</math>; |
− | Сложность вычисления матрицы, обратной произвольной матрице ковариации <math>\Sigma</math> методом Гаусса - <math>O( | + | Сложность вычисления матрицы, обратной произвольной матрице ковариации <math>\Sigma</math> методом Гаусса - <math>O(M^3)</math>; |
− | Сложность предварительного вычисления определителей и обратных матриц матрицам ковариации - <math>O( | + | Сложность предварительного вычисления определителей и обратных матриц матрицам ковариации - <math>O(K * M^3)</math> |
− | С учетом заранее вычисленных определителей и обратных матриц сложность вычисления <math>p(x | \mu,\Sigma)</math> составит: <math>O( | + | С учетом заранее вычисленных определителей и обратных матриц сложность вычисления <math>p(x | \mu,\Sigma)</math> составит: <math>O(M^2)</math>; |
− | Суммарно сложность шага E составит <math>O( | + | Суммарно сложность шага E составит <math>O(N * M^2 * K + K * M^3)</math> |
'''Шаг M''' | '''Шаг M''' | ||
− | Сложность обновления весов кластеров <math>O( | + | Сложность обновления весов кластеров <math>O(K * N)</math>; |
+ | |||
+ | Сложность обновления центров кластеров <math>O(K * M * N)</math>; | ||
+ | |||
+ | Сложность обновления одной матрицы ковариации <math>O(N * M)</math>; | ||
+ | |||
+ | Сложность обновления всех матриц ковариации <math>O(K * M * N)</math>; | ||
+ | |||
+ | Суммарно сложность шага M алгоритма - <math>O(K * M * N)</math>; | ||
+ | |||
+ | |||
+ | |||
+ | Таким образом, сложность одной итерации EM алгоритма для псевдокода, указанного в предыдущем разделе, составит <math>O(K * M^2 * N + K * M^3 + K * M * N)</math>. | ||
− | + | В предположении, что размерность объектов много меньше их количества <math>M << N</math> получим сложность итерации <math>O(K * M^2 * N)</math>. | |
− | + | Нетрудно заметить, что вычисление определителя и обратной матрицы к произвольной матрице ковариации составляет большу́ю часть работы шага E и существенно влияет на его сложность. Поэтому эти вычисления проведены заранее - вынесены за основной цикл шага E<ref name=ParallelEMComparissonGraphLab>Henggang Cui, Jinliang Wei, Wei Dai. Parallel Implementation of Expectation-Maximization for Fast Convergence - School of Computer Science, Carnegie Mellon University. [https://users.ece.cmu.edu/~hengganc/archive/report/final.pdf]</ref>. Помимо этого сложность этих вычислений можно сократить еще сильнее, если ввести ограничения структуру матрицы ковариации <math>\Sigma</math>. | |
− | + | Наиболее часто используемым ограничением, позволяющим эффективно распараллеливать алгоритм<ref>Вячеслав Орешков. EМ — масштабируемый алгоритм кластеризации [https://basegroup.ru/community/articles/em]</ref>, является предположение, что матрица ковариации имеет диагональную форму. Тогда вычисление её определителя и обратной матрицы имеет сложнось <math>O(M)</math>, а суммарная сложность шага E и сложность всей итерации составит <math>O(K * M * N)</math>. | |
− | + | == Информационный граф == | |
+ | Макроструктура алгоритма представлена на рисунке 1 (см. [[#Макроструктура алгоритма|раздел 1.4]] статьи). | ||
+ | На рисунке 2 представлена схема вычислений E-шага алгоритма для одного объекта: для объекта <math>x_{i}</math> и параметров <math>w_{j}, \mu_{j}, \Sigma_{j}, j=1..K</math> по формулам [[#Математическое описание алгоритма|раздела 1.2]] вычисляются ненормированные значения скрытых переменных <math>g_{j}</math>, нормировочный коэффициент ''normСoeff'', и получаются итоговые значения <math>\gamma_{i,j}</math>. | ||
+ | [[file:EM_Expectation_step.png|thumb|center|400px|Рисунок 2: E-шаг EM алгоритма для одного объекта]] | ||
− | + | На рисунке 3 представлена схема вычислений E-шага алгоритма для всех объектов: для объекта <math>x_{j}</math> и параметров-векторов <math>w, \mu, \Sigma</math> по формулам [[#Математическое описание алгоритма|раздела 1.2]] вычисляются значения скрытых переменных. <math>E_{i}, i=1..N</math> - маркооперации, E-шаги для соответствующих объектов. | |
− | + | [[file:EM_Expectation_step_all.png|thumb|center|400px|Рисунок 3: E-шаг EM алгоритма для всех объектов]] | |
− | + | На рисунке 4 представлена схема вычислений M-шага алгоритма для одного кластера: для объектов <math>x_{i}</math> и скрытых переменных <math>\gamma_{i,j}</math>, соответствующих степени принадлежности этих объектов кластеру <math>j</math>, согласно формулам шага M [[#Математическое описание алгоритма|раздела 1.2]] вычисляется нормировочный коэффициент <math>N_{j}</math> и параметры этого кластера: вес <math>w_{j}^{new}</math>, матрица ковариации <math>\Sigma_{j}^{new}</math> и центр кластера <math>\mu_{j}^{new}</math>. | |
− | + | [[file:EM_Maximization_step.png|thumb|center|400px|Рисунок 4: M-шаг EM алгоритма для одного кластера]] | |
− | = | + | На рисунке 5 представлена схема вычислений M-шага алгоритма для всех кластеров. <math>M_{j}, j=1..K</math> соответствует макрооперации - M-шагу для соответствующего кластера и скрытых переменных. |
− | + | ||
− | + | [[file:EM_Maximization_step_all.png|thumb|center|400px|Рисунок 5: M-шаг EM алгоритма для всех кластеров]] | |
+ | |||
+ | == Ресурс параллелизма алгоритма == | ||
+ | Из схем вычислений шагов E и M шагов алгоритма, приведенным в [[#Информационный граф|предыдущем разделе]], видно, что если размерность матрицы ковариации небольшая (тогда вычисление определителя и обратной матрицы матрице ковариации не является трудоемкой операцией), обе макрооперации EM алгоритма можно распараллеливать как по кластерам, так и по объектам. | ||
+ | |||
+ | Как правило количество кластеров, на которые нужно разбить объекты, несравнимо меньше количества объектов, поэтому распараллеливание в большинстве существующих реализаций происходит по объектам: обработка объектов по возможности равномерно распределяются по кластерам. Теоретически (если считать передачу данных между процессорами мгновенной) при наличии неограниченного числа процессоров сложность алгоритма составит <math>O(M^2 * K)</math>. | ||
+ | |||
+ | == Входные и выходные данные алгоритма == | ||
'''Входные данные''': | '''Входные данные''': | ||
− | * Целое неотрицательное число <math> | + | * Целое неотрицательное число <math>K</math> - количество кластеров; |
− | * Значения координат объектов <math>x_{i}</math> - <math> | + | * Значения координат объектов <math>x_{i}</math> - <math>N</math> объектов с <math>M</math> координатами каждый. |
'''Объем входных данных''': | '''Объем входных данных''': | ||
− | * <math> | + | * <math>N * M</math> вещественных чисел (если координаты объектов - вещественные числа), <math>1</math> целое неотрицательное число. |
'''Выходные данные''': | '''Выходные данные''': | ||
− | * Вектор длины <math> | + | * Вектор длины <math>N</math> - для каждого объекта указан номер кластера, к которому он отнесен. Объект относится по результатам работы алгоритма тому кластеру, вероятность принадлежности которому (значение соответствующей скрытой переменной) максимальна. |
'''Объем выходных данных''': | '''Объем выходных данных''': | ||
− | * <math> | + | * <math>N</math> целых неотрицательных чисел. |
− | + | ||
− | EM-алгоритм обладает следующими преимуществами и недостатками: | + | == Свойства алгоритма == |
+ | === Преимущества и недостатки === | ||
+ | EM-алгоритм обладает следующими преимуществами и недостатками<ref>EM алгоритм, его модификации и обобщения. machinelearning.ru [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]</ref>: | ||
'''Преимущества:''' | '''Преимущества:''' | ||
− | + | # Слабо чувствителен к выбросам | |
+ | # Прост в реализации | ||
+ | # Быстро сходится при удачном выборе начальных значений параметров. | ||
+ | |||
+ | '''Недостатки:''' | ||
+ | |||
+ | # Неустойчив к выбору начальных значений параметров - от них зависит как скорость сходимости, так и результат работы: алгоритм может сойтись к локальному экстремуму функции правдоподобия, который может быть существенно ниже глобального. | ||
+ | # Алгоритм не находит оптимальное количество кластеров. Количество кластеров, на которые нужно разбить множество объектов, является параметром алгоритма. | ||
+ | # При большой размерности пространства объектов выдвинутое предположение о модели их распределения может быть некорректным. | ||
− | + | === Влияние структуры матрицы ковариации на форму кластера === | |
− | + | [[file:DifferentCovariance.png|thumb|center|800px|Рисунок 6: Влияние структуры матрицы ковариации на форму кластеров для двумерного случая. a) - матрицы ковариации являются диагональными, b) - матрицы ковариации произвольны. <ref>Mathworks: Clustering Using Gaussian Mixture Models [http://www.mathworks.com/help/stats/clustering-using-gaussian-mixture-models.html]</ref>]] | |
− | '' | + | Рисунок 6 наглядно показывает, что при диагональной форме матриц ковариации полуоси эллипса (в общем случае - эллипсоида), соответствующего распределению объектов, отнесенных кластеру, параллельны осям координат объектов. В противном случае полуоси могут располагаться произвольным образом. |
+ | |||
+ | === Общие свойства === | ||
+ | [[глоссарий#Вычислительная мощность|''Вычислительная мощность'']] EM алгоритма равна <math>M * K * it</math>, где <math>M</math> - размерность объекта, <math>K</math> - количество кластеров, а <math>it</math> - количество итераций. | ||
− | + | EM-алгоритм является итерационным, и число итераций в общем случае не фиксировано: оно зависит от условий остановки, начальных приближений и расположения объектов. Поэтому этот алгоритм ''недетерминированный''. По тем же причинам EM-алгоритм является ''неустойчивым''. | |
− | + | = Программная реализация алгоритма = | |
− | 3. | + | == Особенности реализации последовательного алгоритма == |
+ | === Инициализация параметров распределений === | ||
+ | Наиболее распространены 3 способа инициализации центров распределений для EM алгоритма: | ||
+ | # центр распределения совпадает со случайно выбранным объектом <math>x_{i}</math> | ||
+ | # все наблюдения распределяются по кластерам случайным образом, а затем для каждого кластера пересчитывается цент соответствующего распределения | ||
+ | # центры кластеров получаются в результате работы алгоритма k-means. Использование такой инициализации оправдано тем, что k-means требует меньших вычислительных затрат и быстрее сходится<ref>EM алгоритм (пример) machinelearning.ru [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])</ref> | ||
− | == | + | == Локальность данных и вычислений == |
− | == | + | == Возможные способы и особенности параллельной реализации алгоритма == |
− | + | == Масштабируемость алгоритма и его реализации == | |
− | |||
− | |||
− | |||
− | |||
− | + | Многие статьи<ref name=FREERIDE>Leonid Glimcher, Gagan Agrawal. Parallelising EM Algorithm on a Cluster of SMPs [http://web.cse.ohio-state.edu/~agrawal/allpapers/europar04.pdf]</ref><ref name=CV>López-de-Teruel P. E., García J. M., Acacio; The Parallel EM Algorithm and its Applications in Computer Vision [http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.4.8470&rep=rep1&type=pdf]</ref><ref name=MULTITHREAD>Sharon X. Lee, Kaleb L. Leemaqz, Geoffrey J. McLachlan. A simple multithreaded implementation of the EM algorithm for mixture models [https://arxiv.org/pdf/1606.02054.pdf]</ref><ref name=CUDA>Araújo, G.F., H.T. Macedo, M.T. Chella, C.A.E. Montesco and M.V.O. Medeiros. PARALLEL IMPLEMENTATION OF EXPECTATION-MAXIMISATION ALGORITHM FOR THE TRAINING OF GAUSSIAN MIXTURE MODELS [https://ri.ufs.br/bitstream/123456789/1764/1/ExpectationMaximisationAlgorithm.pdf]</ref>, в которых представлены параллельные реализации EM алгоритма, описывают помимо использующихся программных и аппаратных решений результаты проведенных экспериментов с данными. В качестве примера в этом разделе будут рассмотрены результаты эксперементов, проведенных со связующим программным обеспечением (''middleware'') FREERIDE<ref>Leonid Glimchera, Ruoming Jinb, Gagan Agrawala. Middleware for data mining applications on clusters and grids [http://www.cs.kent.edu/~jin/Papers/JPDC08.pdf]</ref> (FRamework for Rapid Implementation of Datamining Engines), обеспечивающим возможность распараллеливания как для систем с общей памятью (в частности, симметричных мультипроцессорных систем), так и для систем с распределенной памятью. | |
− | + | Эксперименты проводились на кластере машин Pentium 700MHz связанных при помощи Myrinet<ref>https://en.wikipedia.org/wiki/Myrinet</ref> LANai 7.0 на трех наборах данных: | |
− | + | '''а)''' 3.3 * 10^6 объектов | |
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | + | '''b)''' 6.2 * 10^6 объектов | |
− | + | '''c)''' 12.5 * 10^6 объектов | |
− | + | Объекты представляют собой координаты точек 10-мерного пространства. | |
+ | |||
+ | Количество машин в кластере варьировалось от 1 до 8, при этом на каждой машине запускалось до четырех параллельных потоков. | ||
− | = | + | Результаты экспериментов представлены на рисунках 7-9<ref name=FREERIDE/>: |
− | == | + | {|align="center" |
+ | |-valign="top" | ||
+ | |[[file:EM_FREERIDE_experiments_a.png|thumb|center|350px|Рисунок 7: набор данных a)]] | ||
+ | |[[file:EM_FREERIDE_experiments_b.png|thumb|center|357px|Рисунок 8: набор данных b)]] | ||
+ | |[[file:EM_FREERIDE_experiments_c.png|thumb|center|351px|Рисунок 9: набор данных c)]] | ||
+ | |} | ||
− | + | Несмотря на то, что используемые в данном исследовании машины, как и коммуникационная среда, сейчас не используются и считаются устаревшими, результат приведенных экспериментов наглядно показывает, что EM алгоритм хорошо параллелится, зависимость ускорения от числа узлов в кластере с небольшим количеством узлов близка к линейной, а также, что скорость работы алгоритма линейно зависит от числа объектов. | |
− | + | '''''Данный раздел может быть дополнен''''' | |
− | == | + | == Динамические характеристики и эффективность реализации алгоритма == |
− | |||
− | + | == Выводы для классов архитектур == | |
− | + | == Существующие реализации алгоритма == | |
+ | Примеры библиотек, содержащих последовательную реализацию EM-алгоритма, или отдельных реализаций: | ||
− | [ | + | # [http://docs.opencv.org/2.4/modules/ml/doc/expectation_maximization.html OpenCV(C++)] ([https://github.com/opencv/opencv/blob/master/modules/ml/src/em.cpp github]) |
+ | # [http://www.mathworks.com/help/stats/clustering-using-gaussian-mixture-models.html Matlab] | ||
+ | # [http://scikit-learn.org/stable/modules/mixture.html Scikit learn (Python)] | ||
+ | # [http://mlpack.org/index.html MLPack(C++)] | ||
− | + | Примеры библиотек, содержащих параллельную реализацию EM-алгоритма, или отдельных реализаций: | |
− | + | ''Здесь приводятся как конкретные реализации, так и статьи, описывающие подходы к реализациям распараллеливания EM-алгоритма'' | |
− | [ | + | # [https://spark.apache.org/ Apache Spark] |
+ | # [http://accord-framework.net/ Accord.NET Framework (C#)] | ||
+ | # Статья об использовании FREERIDE middleware для распараллеливания EM алгоритма <ref name=FREERIDE/> | ||
+ | # [http://andrewharp.com/gmmcuda реализация с использованием CUDA] | ||
+ | # [https://github.com/Corv/CUDA-GMM-MPI еще одна реализация с использованием CUDA] | ||
+ | # статья об использовании параллельного EM алгоритма в задачах компьютерного зрения <ref name=CV/> | ||
+ | # статья о распараллеливании EM алгоритма с использованием CUDA <ref name=CUDA/> | ||
− | + | = Литература = | |
+ | <references /> |
Текущая версия на 16:16, 21 ноября 2016
EM-алгоритм кластеризации | |
Последовательный алгоритм | |
Последовательная сложность | [math]O(iter * N * M^2 * K)[/math] |
Объём входных данных | [math]N * M + 1[/math] |
Объём выходных данных | [math]N[/math] |
Параллельный алгоритм | |
Высота ярусно-параллельной формы | [math]O(iter * M^2 * K)[/math] |
Ширина ярусно-параллельной формы | [math]O(iter * N)[/math] |
Автор описания: Зинченко Д.А.
Содержание
- 1 Свойства и структура алгоритма
- 1.1 Общее описание алгоритма
- 1.2 Математическое описание алгоритма
- 1.3 Вычислительное ядро алгоритма
- 1.4 Макроструктура алгоритма
- 1.5 Схема реализации последовательного алгоритма
- 1.6 Последовательная сложность алгоритма
- 1.7 Информационный граф
- 1.8 Ресурс параллелизма алгоритма
- 1.9 Входные и выходные данные алгоритма
- 1.10 Свойства алгоритма
- 2 Программная реализация алгоритма
- 2.1 Особенности реализации последовательного алгоритма
- 2.2 Локальность данных и вычислений
- 2.3 Возможные способы и особенности параллельной реализации алгоритма
- 2.4 Масштабируемость алгоритма и его реализации
- 2.5 Динамические характеристики и эффективность реализации алгоритма
- 2.6 Выводы для классов архитектур
- 2.7 Существующие реализации алгоритма
- 3 Литература
1 Свойства и структура алгоритма
1.1 Общее описание алгоритма
В задачах оптимизации EM-алгоритмом называют итеративную процедуру поиска численного значения экстремума какой-либо функции. В частности, в статистике этот алгоритм используется для оценки максимального правдоподобия. Впервые название EM-алгоритм было предложено в 1977 году[1], однако его идеи были описаны и раньше[2].
EM-алгоритм применяется для решения задач двух типов [2]:
- Анализ неполных данных (данных с пропусками).
- Оценка максимального правдоподобия в случае, если функцию правдоподобия трудно исследовать аналитически. В этом случае введение набора скрытых переменных может существенно упростить задачу.
Задача кластеризации относится к задачам второго типа. В данной статье рассматривается следующая постановка задачи кластеризации:
Дано множество из N объектов, каждый из которых представляет собой точку в M-мерном метрическом пространстве. Каждому измерению пространства соответствует некоторое свойство объекта. Необходимо разбить это множество на K подмножеств (кластеров) так, чтобы, чтобы элементы каждого подмножества существенно отличались по свойствам от элементов других подмножеств.
Использование алгоритма EM подразумевает следующее предположение об объектах:
Пусть объекты, которые необходимо разбить на кластеры, появляются случайным образом и независимо друг от друга согласно вероятностному распределению, равному смеси (линейной комбинации) распределений кластеров. В данной статье в дальнейшем будет считаться, что распределение каждого кластера является многомерным нормальным распределением с произвольной матрицей ковариации - наиболее часто используемое предположение.
Тогда каждый объект принадлежит каждому кластеру, но с разной вероятностью. EM-алгоритм итеративно оценивает параметры распределений кластеров, максимизируя логарифмическую функцию максимального правдоподобия. После окончания работы алгоритма объект будет отнесен к кластеру, вероятность принадлежности которому максимальна. Вводимый набор скрытых переменных для каждого объекта - вероятности того, что объект принадлежит каждому из кластеров.
Итерация алгоритма состоит из двух последовательных шагов.
- (Expectation) Вычисляются новые ожидаемые значения скрытых переменных.
- (Maximization) Решается задача максимизации правдоподобия: По текущим значениям скрытых переменных обновляются параметры распределений для кластеров: математическое ожидание, дисперсия и вероятность появления объектов из кластера.
1.2 Математическое описание алгоритма
Используемые обозначения:
[math]N[/math] - количество объектов, [math]K[/math] - количество кластеров, [math]M[/math] - количество параметров объекта
[math]\gamma_{n,k}[/math] - значения скрытых переменных
[math]N_{k}[/math] - количество элементов отнесенных к кластеру [math]k[/math]
[math]x_{1},...,x_N[/math] - объекты
Параметры распределений:
[math]w_{1},...,w_{K}[/math] - вероятности появления объектов из кластеров (веса кластеров), [math]\sum_{k=1}^Kw_{k}=1[/math]
[math]\mu_{1},...,\mu_{K}[/math] - центры гауссиан кластеров
[math]\Sigma_{1},...,\Sigma_{K}[/math] - матрицы ковариации гауссиан кластеров, каждая матрица имеет размерность [math]M[/math]x[math]M[/math]
Смесь распределений:
Цель алгоритма - вычислить параметры распределений, максимизирующих логарифм функции правдоподобия:
- В начале работы алгоритма задаются параметры начального приближения:
Поскольку от начального приближения может сильно зависеть результат работы алгоритма, вместо инициализации кластеров случайными центрами, равными весами и диагональными матрицами ковариации, как в приведенных формулах, часто применяют другие способы (подробнее - в особенностях реализации)
- Далее итеративно выполняется следующая пара процедур:
- E-шаг: используя текущее значение параметров [math]w_{1},...,w_{k};\mu_{1},...,\mu_{k};\Sigma^{1},...,\Sigma^{k}[/math], вычисляем значение вектора скрытых переменных [math]\gamma[/math]:
[math]\gamma_{n,k}=\frac{w_{k}p(x_{n}|\mu_{k},\Sigma^{k})}{\sum_{j=1}^K w_{j}p(x_{n}|\mu_{j},\Sigma^{j})}, [/math] знаменатель дроби - нормировочный коэффициент [math]p(x|\mu,\Sigma)=\frac{1}{{(2\pi)}^{\frac{N}{2}}\sqrt{|\Sigma|}}exp\biggl\{-\frac{1}{2}{(x-\mu)}^T{\Sigma}^{-1}(x-\mu)\biggr\}[/math] - плотность [math]N[/math]-мерного нормального распределения - М-шаг: переоценка параметров по текущим значениям скрытых переменных для [math]k = 1..K[/math]:
- E-шаг: используя текущее значение параметров [math]w_{1},...,w_{k};\mu_{1},...,\mu_{k};\Sigma^{1},...,\Sigma^{k}[/math], вычисляем значение вектора скрытых переменных [math]\gamma[/math]:
В EM алгоритме используется один из следующих критериев остановки:
- Норма разности векторов скрытых переменных на текущей итерации не превышает заданную константу
- Прошествие заданного количества итераций
- Изменение логарифмического правдоподобия меньше заданной константы.
1.3 Вычислительное ядро алгоритма
Вычислительное ядро алгоритма состоит из двух шагов, повторяющихся каждую итерацию:
- Обновление скрытых переменных на основе текущих параметров распределений
- Обновление параметров кластеров: весов, центров и матриц ковариации
1.4 Макроструктура алгоритма
Как записано в описании вычислительного ядра алгоритма, основную часть алгоритма составляют итеративно проводимые макрооперации: Expectation шаг - вычисление новых значений скрытых переменных, Maximization шаг - обновление параметров кластеров. Макроструктура алгоритма представлена на рисунке 1:
1.5 Схема реализации последовательного алгоритма
Примерная схема реализации EM алгоритма с предварительной обработкой матриц ковариации:
Инициализировать веса и центры кластеров;
Инициализировать матрицы ковариации;
// Основной цикл
Пока условие остановки не выполнено
// E-шаг
В цикле для каждого кластера
Вычислить определитель матрицы ковариации Sigma;
Вычислить обратную матрицу матрице Sigma;
В цикле для каждого объекта
В цикле для каждого кластера
Вычислить и запомнить p(x | mu, Sigma) для текущего объекта и кластера;
Просуммировать вычисленные p(x | mu, Sigma);
В цикле для каждого кластера
Обновить скрытые переменные gamma;
// M-шаг
Обновить веса w;
Обновить центры кластеров mu;
Обновить матрицы ковариации Sigma;
1.6 Последовательная сложность алгоритма
Количество итераций алгоритма зависит от условия остановки, поэтому будем рассматривать сложность одной итерации.
Шаг E
Сложность вычисления определителя произвольной матрицы ковариации [math]\Sigma[/math] методом Гаусса - [math]O(M^3)[/math];
Сложность вычисления матрицы, обратной произвольной матрице ковариации [math]\Sigma[/math] методом Гаусса - [math]O(M^3)[/math];
Сложность предварительного вычисления определителей и обратных матриц матрицам ковариации - [math]O(K * M^3)[/math]
С учетом заранее вычисленных определителей и обратных матриц сложность вычисления [math]p(x | \mu,\Sigma)[/math] составит: [math]O(M^2)[/math];
Суммарно сложность шага E составит [math]O(N * M^2 * K + K * M^3)[/math]
Шаг M
Сложность обновления весов кластеров [math]O(K * N)[/math];
Сложность обновления центров кластеров [math]O(K * M * N)[/math];
Сложность обновления одной матрицы ковариации [math]O(N * M)[/math];
Сложность обновления всех матриц ковариации [math]O(K * M * N)[/math];
Суммарно сложность шага M алгоритма - [math]O(K * M * N)[/math];
Таким образом, сложность одной итерации EM алгоритма для псевдокода, указанного в предыдущем разделе, составит [math]O(K * M^2 * N + K * M^3 + K * M * N)[/math].
В предположении, что размерность объектов много меньше их количества [math]M \lt \lt N[/math] получим сложность итерации [math]O(K * M^2 * N)[/math].
Нетрудно заметить, что вычисление определителя и обратной матрицы к произвольной матрице ковариации составляет большу́ю часть работы шага E и существенно влияет на его сложность. Поэтому эти вычисления проведены заранее - вынесены за основной цикл шага E[3]. Помимо этого сложность этих вычислений можно сократить еще сильнее, если ввести ограничения структуру матрицы ковариации [math]\Sigma[/math].
Наиболее часто используемым ограничением, позволяющим эффективно распараллеливать алгоритм[4], является предположение, что матрица ковариации имеет диагональную форму. Тогда вычисление её определителя и обратной матрицы имеет сложнось [math]O(M)[/math], а суммарная сложность шага E и сложность всей итерации составит [math]O(K * M * N)[/math].
1.7 Информационный граф
Макроструктура алгоритма представлена на рисунке 1 (см. раздел 1.4 статьи).
На рисунке 2 представлена схема вычислений E-шага алгоритма для одного объекта: для объекта [math]x_{i}[/math] и параметров [math]w_{j}, \mu_{j}, \Sigma_{j}, j=1..K[/math] по формулам раздела 1.2 вычисляются ненормированные значения скрытых переменных [math]g_{j}[/math], нормировочный коэффициент normСoeff, и получаются итоговые значения [math]\gamma_{i,j}[/math].
На рисунке 3 представлена схема вычислений E-шага алгоритма для всех объектов: для объекта [math]x_{j}[/math] и параметров-векторов [math]w, \mu, \Sigma[/math] по формулам раздела 1.2 вычисляются значения скрытых переменных. [math]E_{i}, i=1..N[/math] - маркооперации, E-шаги для соответствующих объектов.
На рисунке 4 представлена схема вычислений M-шага алгоритма для одного кластера: для объектов [math]x_{i}[/math] и скрытых переменных [math]\gamma_{i,j}[/math], соответствующих степени принадлежности этих объектов кластеру [math]j[/math], согласно формулам шага M раздела 1.2 вычисляется нормировочный коэффициент [math]N_{j}[/math] и параметры этого кластера: вес [math]w_{j}^{new}[/math], матрица ковариации [math]\Sigma_{j}^{new}[/math] и центр кластера [math]\mu_{j}^{new}[/math].
На рисунке 5 представлена схема вычислений M-шага алгоритма для всех кластеров. [math]M_{j}, j=1..K[/math] соответствует макрооперации - M-шагу для соответствующего кластера и скрытых переменных.
1.8 Ресурс параллелизма алгоритма
Из схем вычислений шагов E и M шагов алгоритма, приведенным в предыдущем разделе, видно, что если размерность матрицы ковариации небольшая (тогда вычисление определителя и обратной матрицы матрице ковариации не является трудоемкой операцией), обе макрооперации EM алгоритма можно распараллеливать как по кластерам, так и по объектам.
Как правило количество кластеров, на которые нужно разбить объекты, несравнимо меньше количества объектов, поэтому распараллеливание в большинстве существующих реализаций происходит по объектам: обработка объектов по возможности равномерно распределяются по кластерам. Теоретически (если считать передачу данных между процессорами мгновенной) при наличии неограниченного числа процессоров сложность алгоритма составит [math]O(M^2 * K)[/math].
1.9 Входные и выходные данные алгоритма
Входные данные:
- Целое неотрицательное число [math]K[/math] - количество кластеров;
- Значения координат объектов [math]x_{i}[/math] - [math]N[/math] объектов с [math]M[/math] координатами каждый.
Объем входных данных:
- [math]N * M[/math] вещественных чисел (если координаты объектов - вещественные числа), [math]1[/math] целое неотрицательное число.
Выходные данные:
- Вектор длины [math]N[/math] - для каждого объекта указан номер кластера, к которому он отнесен. Объект относится по результатам работы алгоритма тому кластеру, вероятность принадлежности которому (значение соответствующей скрытой переменной) максимальна.
Объем выходных данных:
- [math]N[/math] целых неотрицательных чисел.
1.10 Свойства алгоритма
1.10.1 Преимущества и недостатки
EM-алгоритм обладает следующими преимуществами и недостатками[5]:
Преимущества:
- Слабо чувствителен к выбросам
- Прост в реализации
- Быстро сходится при удачном выборе начальных значений параметров.
Недостатки:
- Неустойчив к выбору начальных значений параметров - от них зависит как скорость сходимости, так и результат работы: алгоритм может сойтись к локальному экстремуму функции правдоподобия, который может быть существенно ниже глобального.
- Алгоритм не находит оптимальное количество кластеров. Количество кластеров, на которые нужно разбить множество объектов, является параметром алгоритма.
- При большой размерности пространства объектов выдвинутое предположение о модели их распределения может быть некорректным.
1.10.2 Влияние структуры матрицы ковариации на форму кластера
Рисунок 6 наглядно показывает, что при диагональной форме матриц ковариации полуоси эллипса (в общем случае - эллипсоида), соответствующего распределению объектов, отнесенных кластеру, параллельны осям координат объектов. В противном случае полуоси могут располагаться произвольным образом.
1.10.3 Общие свойства
Вычислительная мощность EM алгоритма равна [math]M * K * it[/math], где [math]M[/math] - размерность объекта, [math]K[/math] - количество кластеров, а [math]it[/math] - количество итераций.
EM-алгоритм является итерационным, и число итераций в общем случае не фиксировано: оно зависит от условий остановки, начальных приближений и расположения объектов. Поэтому этот алгоритм недетерминированный. По тем же причинам EM-алгоритм является неустойчивым.
2 Программная реализация алгоритма
2.1 Особенности реализации последовательного алгоритма
2.1.1 Инициализация параметров распределений
Наиболее распространены 3 способа инициализации центров распределений для EM алгоритма:
- центр распределения совпадает со случайно выбранным объектом [math]x_{i}[/math]
- все наблюдения распределяются по кластерам случайным образом, а затем для каждого кластера пересчитывается цент соответствующего распределения
- центры кластеров получаются в результате работы алгоритма k-means. Использование такой инициализации оправдано тем, что k-means требует меньших вычислительных затрат и быстрее сходится[7]
2.2 Локальность данных и вычислений
2.3 Возможные способы и особенности параллельной реализации алгоритма
2.4 Масштабируемость алгоритма и его реализации
Многие статьи[8][9][10][11], в которых представлены параллельные реализации EM алгоритма, описывают помимо использующихся программных и аппаратных решений результаты проведенных экспериментов с данными. В качестве примера в этом разделе будут рассмотрены результаты эксперементов, проведенных со связующим программным обеспечением (middleware) FREERIDE[12] (FRamework for Rapid Implementation of Datamining Engines), обеспечивающим возможность распараллеливания как для систем с общей памятью (в частности, симметричных мультипроцессорных систем), так и для систем с распределенной памятью.
Эксперименты проводились на кластере машин Pentium 700MHz связанных при помощи Myrinet[13] LANai 7.0 на трех наборах данных:
а) 3.3 * 10^6 объектов
b) 6.2 * 10^6 объектов
c) 12.5 * 10^6 объектов
Объекты представляют собой координаты точек 10-мерного пространства.
Количество машин в кластере варьировалось от 1 до 8, при этом на каждой машине запускалось до четырех параллельных потоков.
Результаты экспериментов представлены на рисунках 7-9[8]:
Несмотря на то, что используемые в данном исследовании машины, как и коммуникационная среда, сейчас не используются и считаются устаревшими, результат приведенных экспериментов наглядно показывает, что EM алгоритм хорошо параллелится, зависимость ускорения от числа узлов в кластере с небольшим количеством узлов близка к линейной, а также, что скорость работы алгоритма линейно зависит от числа объектов.
Данный раздел может быть дополнен
2.5 Динамические характеристики и эффективность реализации алгоритма
2.6 Выводы для классов архитектур
2.7 Существующие реализации алгоритма
Примеры библиотек, содержащих последовательную реализацию EM-алгоритма, или отдельных реализаций:
Примеры библиотек, содержащих параллельную реализацию EM-алгоритма, или отдельных реализаций:
Здесь приводятся как конкретные реализации, так и статьи, описывающие подходы к реализациям распараллеливания EM-алгоритма
- Apache Spark
- Accord.NET Framework (C#)
- Статья об использовании FREERIDE middleware для распараллеливания EM алгоритма [8]
- реализация с использованием CUDA
- еще одна реализация с использованием CUDA
- статья об использовании параллельного EM алгоритма в задачах компьютерного зрения [9]
- статья о распараллеливании EM алгоритма с использованием CUDA [11]
3 Литература
- ↑ 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.
- ↑ 2,0 2,1 В.Ю.Королёв. ЕМ-алгоритм, его модификации и их применение к задаче разделения смесей вероятностных распределений.
- ↑ Henggang Cui, Jinliang Wei, Wei Dai. Parallel Implementation of Expectation-Maximization for Fast Convergence - School of Computer Science, Carnegie Mellon University. [1]
- ↑ Вячеслав Орешков. EМ — масштабируемый алгоритм кластеризации [2]
- ↑ EM алгоритм, его модификации и обобщения. machinelearning.ru [3]
- ↑ Mathworks: Clustering Using Gaussian Mixture Models [4]
- ↑ EM алгоритм (пример) machinelearning.ru [5])
- ↑ 8,0 8,1 8,2 Leonid Glimcher, Gagan Agrawal. Parallelising EM Algorithm on a Cluster of SMPs [6]
- ↑ 9,0 9,1 López-de-Teruel P. E., García J. M., Acacio; The Parallel EM Algorithm and its Applications in Computer Vision [7]
- ↑ Sharon X. Lee, Kaleb L. Leemaqz, Geoffrey J. McLachlan. A simple multithreaded implementation of the EM algorithm for mixture models [8]
- ↑ 11,0 11,1 Araújo, G.F., H.T. Macedo, M.T. Chella, C.A.E. Montesco and M.V.O. Medeiros. PARALLEL IMPLEMENTATION OF EXPECTATION-MAXIMISATION ALGORITHM FOR THE TRAINING OF GAUSSIAN MIXTURE MODELS [9]
- ↑ Leonid Glimchera, Ruoming Jinb, Gagan Agrawala. Middleware for data mining applications on clusters and grids [10]
- ↑ https://en.wikipedia.org/wiki/Myrinet