Участник:Noite/EM-алгоритм кластеризации
Эта работа успешно выполнена Преподавателю: в основное пространство, в подстраницу Данное задание было проверено и зачтено. Проверено Dexter и Преподаватель 1. |
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