Участник:Александр Куваев/Алгоритм кластеризации, основанный на максимизации ожидания: различия между версиями
ASA (обсуждение | вклад) |
|||
(не показано 46 промежуточных версий 3 участников) | |||
Строка 1: | Строка 1: | ||
− | Авторы описания: Куваев А.С. | + | {{Assignment|Algoman|Dexter}} |
+ | |||
+ | Авторы описания: | ||
+ | * Куваев А.С., группа 620 (математическая постановка задачи, схема реализации алгоритма, оценка сложности, описание ресурса параллелизма, реализация алгоритма) | ||
+ | * Щенявская Е.В., группа 616 (описание алгоритма, визуализация информационного графа и примера работы алгоритма, описание свойств алгоритма, замеры времени выполнения и визуализация результата) | ||
== Свойства и структура алгоритмов == | == Свойства и структура алгоритмов == | ||
=== Общее описание алгоритма === | === Общее описание алгоритма === | ||
− | Задача кластеризации заключается в разбиении входного множества объектов на непересекающиеся подмножества, называемые кластерами, так, чтобы каждый кластер состоял из схожих объектов, а объекты разных кластеров существенно отличались друг от друга. | + | Задача кластеризации заключается в разбиении входного множества объектов на непересекающиеся подмножества, называемые кластерами, так, чтобы каждый кластер состоял из схожих объектов, а объекты разных кластеров существенно отличались друг от друга<ref name=vorontsov>Воронцов К.В., Математические методы обучения по прецедентам.</ref>. |
− | Решение этой задачи принципиально неоднозначно по следующим причинам: | + | Решение этой задачи принципиально неоднозначно по следующим причинам<ref name=vorontsov/>: |
# результат кластеризации зависит от способа задания меры сходства объектов выборки | # результат кластеризации зависит от способа задания меры сходства объектов выборки | ||
# не существует однозначно наилучшего критерия качества кластеризации | # не существует однозначно наилучшего критерия качества кластеризации | ||
Строка 13: | Строка 17: | ||
По описанной выше причине существует большое число алгоритмов кластеризации, приводящих к различным разбиениям исходного множества объектов. На этой странице представлено описание одного из таких методов — EM-алгоритма. EM-алгоритм опирается на предположение о вероятностной природе данных: элементы выборки получены случайно и независимо из смеси распределений с фиксированным числом компонент <math>k</math>. Таким образом, плотность распределения на множестве объектов имеет следующий вид: | По описанной выше причине существует большое число алгоритмов кластеризации, приводящих к различным разбиениям исходного множества объектов. На этой странице представлено описание одного из таких методов — EM-алгоритма. EM-алгоритм опирается на предположение о вероятностной природе данных: элементы выборки получены случайно и независимо из смеси распределений с фиксированным числом компонент <math>k</math>. Таким образом, плотность распределения на множестве объектов имеет следующий вид: | ||
− | <math>p(x)=\sum_{j=1}^{k}w_{j}p_{j}(x), \sum_{j=1}^{k}w_{j}=1, w_{j} \ge 0 </math>, где <math>p_{j}</math> - плотность распределения <math>j</math>-й компоненты смеси. | + | :<math>p(x)=\sum_{j=1}^{k}w_{j}p_{j}(x), \ \sum_{j=1}^{k}w_{j}=1, \ w_{j} \ge 0 </math>, где <math>p_{j}</math> - плотность распределения <math>j</math>-й компоненты смеси (кластера). |
− | Везде в дальнейшем будем предполагать, что <math>p_{j}</math> имеют вид многомерных нормальных плотностей с произвольной матрицей ковариации: смеси нормальных распределений позволяют аппроксимировать произвольные непрерывные функции плотности с наперед заданной точностью. Результатом работы EM-алгоритма являются оценки априорных вероятностей компонент смеси <math>w_{j}</math>, а также оценки векторов математических ожиданий и матриц | + | Везде в дальнейшем будем предполагать, что <math>p_{j}</math> имеют вид многомерных нормальных плотностей с произвольной матрицей ковариации: смеси нормальных распределений позволяют аппроксимировать произвольные непрерывные функции плотности с наперед заданной точностью<ref name=korolev>Королёв В.Ю., ЕМ-алгоритм, его модификации и их применение к задаче разделения смесей вероятностных распределений. Теоретический обзор. - М.: ИПИРАН, 2007.</ref>. Результатом работы EM-алгоритма являются оценки априорных вероятностей компонент смеси <math>w_{j}</math>, а также оценки векторов математических ожиданий и матриц ковариаций для каждой компоненты. Зная параметры распределения смеси, каждому объекту будет сопоставляться тот кластер, вероятность принадлежности к которому будет максимальной. |
EM-алгоритм позволяет значительно упростить задачу максимизации правдоподобия выборки путем искусственного введения вспомогательной матрицы скрытых переменных <math>G</math>. Алгоритм заключается в последовательном повторении шагов E(expectation) и M(maximization): | EM-алгоритм позволяет значительно упростить задачу максимизации правдоподобия выборки путем искусственного введения вспомогательной матрицы скрытых переменных <math>G</math>. Алгоритм заключается в последовательном повторении шагов E(expectation) и M(maximization): | ||
Строка 24: | Строка 28: | ||
=== Математическое описание алгоритма === | === Математическое описание алгоритма === | ||
− | Пусть заданы <math>l</math> объектов <math>x_{1}, | + | Пусть заданы <math>l</math> объектов <math>x_{1},\dotsc,x_{l}</math>, каждый из которых описывается <math>n</math> числовыми признаками. Таким образом, определена матрица объектов-признаков <math>X \in \R^{l \times n}.</math> Предполагается, что объекты выбраны случайно и независимо из смеси <math>n</math>-мерных нормальных распределений с фиксированным числом компонент <math>k</math>. Плотность распределения на множестве объектов имеет вид: |
+ | |||
+ | :<math>p(x)=\sum_{j=1}^{k}w_{j}p_{j}(x), \ \sum_{j=1}^{k}w_{j}=1, \ w_{j} \ge 0 </math>, где <math>p_{j}</math> - плотность распределения <math>j</math>-й компоненты смеси. | ||
− | <math> | + | Каждая компонента смеси описывается <math>n</math>-мерным вектором математических ожиданий <math>\mu_{j}</math> и матрицей ковариаций <math>\Sigma_{j}</math> порядка <math>n</math>, <math>j = 1,\dotsc,k</math>. |
− | + | '''Входные данные:''' матрица объектов-признаков <math>X</math>, число кластеров <math>k</math>, максимальное число итераций <math>imax</math>, минимальная величина изменения логарифма правдоподобия <math>\varepsilon</math>. | |
− | ''' | + | '''Выходные данные:''' набор оценок весов, математических ожиданий и ковариационных матриц компонент смеси <math>\theta = (w_{1},...,w_{k}; \; \mu_{1},...,\mu_{k}; \; \Sigma_{1},...,\Sigma_{k})</math>, максимизирующий правдоподобие выборки. |
− | ''' | + | '''Структура алгоритма:''' |
+ | * '''Инициализация параметров:''' существует большое число вариантов инициализации параметров распределения. Один из возможных подходов — задание векторов математических ожиданий компонент случайными элементами выборки, задание матриц ковариаций единичными матрицами и задание весов компонент равными <math>\frac{1}{k}.</math> | ||
+ | * '''Последовательное выполнение шагов E и M''' до тех пор, пока правдоподобие выборки не стабилизируется или не будет достигнуто максимальное число итераций: | ||
+ | ** '''Шаг E:''' вычисление значений скрытых переменных по формуле Байеса: | ||
+ | *:<math>g_{ij} = \frac{w_{j} p_{j}(x_{i})}{\sum_{s=1}^{k} w_{s} p_{s}(x_{i})}</math> , где <math>p_{j}(x_{i}) = \frac{1}{(2\pi)^{n/2}\sqrt{|\Sigma_{j}|}} \exp \biggl( -\frac{1}{2}(x_{i} - \mu_{j})^{T} \Sigma_{j}^{-1} (x_{i} - \mu_{j}) \biggr), \ i = 1,\dotsc,l; \ j = 1,\dotsc,k </math> | ||
+ | ** '''Шаг M:''' перерасчет параметров смеси на основе текущего приближения и матрицы скрытых переменных: | ||
+ | *:<math>w_{j} = \frac{1}{l} \sum_{i=1}^{l} g_{ij}, \ j = 1,\dotsc,k;</math> | ||
+ | *:<math>\mu_{j} = \frac{1}{l w_{j}} \sum_{i=1}^{l} g_{ij} x_{i}, \ j = 1,\dotsc,k;</math> | ||
+ | *:<math>\Sigma_{j} = \frac{1}{l w_{j}} \sum_{i=1}^{l} g_{ij}(x_{i} - \mu_{j})(x_{i} - \mu_{j})^T, \ j = 1,\dotsc,k.</math> | ||
=== Вычислительное ядро алгоритма === | === Вычислительное ядро алгоритма === | ||
+ | Вычислительное ядро EM-алгоритма — процедура последовательного выполнения шагов E и M: | ||
+ | * На шаге E на основе текущего приближения параметров смеси вычисляются ожидаемые значения скрытых переменных | ||
+ | * На шаге M вычисляется следующее приближение параметров смеси на основе текущего приближения и матрицы скрытых переменных | ||
+ | Наиболее трудоемкой операцией с вычислительной точки зрения является шаг E, в ходе которого производится обращение ковариационных матриц, вычисление их определителей, а также многократное перемножение векторов и матриц при вычислении скрытых переменных. | ||
=== Макроструктура алгоритма === | === Макроструктура алгоритма === | ||
+ | Алгоритм состоит из двух итерационно проводимых макроопераций: E-шага и M-шага. | ||
+ | |||
+ | В ходе E-шага используется: | ||
+ | * Обращение ковариационной матрицы и вычисление ее определителя | ||
+ | * Умножение вектор-строки на матрицу и вектор-строки на вектор-столбец | ||
+ | В ходе M-шага используется: | ||
+ | * Умножение вектор-столбца на вектор-строку | ||
+ | * Взвешенное суммирование векторов и матриц | ||
=== Схема реализации последовательного алгоритма === | === Схема реализации последовательного алгоритма === | ||
+ | # Для всех <math>j = 1,\dotsc,k</math>: | ||
+ | #: Инициализировать <math>w_{j}, \ \mu_{j}, \ \Sigma_{j}</math> | ||
+ | # Для всех <math>iter = 1,\dotsc,imax</math>: | ||
+ | ## Шаг E: | ||
+ | ##: Для всех <math>j = 1,\dotsc,k</math>: | ||
+ | ##:: Вычислить <math>|\Sigma_{j}|, \ \Sigma_{j}^{-1}</math> | ||
+ | ##: Для всех <math>i = 1,\dotsc,l</math>: | ||
+ | ##:: <math>sum = 0</math> | ||
+ | ##:: Для всех <math>j = 1,\dotsc,k</math>: | ||
+ | ##::: Вычислить <math>p_{j}(x_{i}) = \frac{1}{(2\pi)^{n/2}\sqrt{|\Sigma_{j}|}} \exp \biggl( -\frac{1}{2}(x_{i} - \mu_{j})^{T} \Sigma_{j}^{-1} (x_{i} - \mu_{j}) \biggr)</math> | ||
+ | ##::: <math>sum = sum + w_{j}p_{j}(x_{i})</math> | ||
+ | ##:: Для всех <math>j = 1,\dotsc,k</math>: | ||
+ | ##::: Вычислить <math>g_{ij} = \frac{w_{j} p_{j}(x_{i})}{sum}</math> | ||
+ | ## Шаг M: | ||
+ | ##: Для всех <math>j = 1,\dotsc,k</math>: | ||
+ | ##:: Вычислить <math>w_{j} = \frac{1}{l} \sum_{i=1}^{l} g_{ij}</math> | ||
+ | ##:: Вычислить <math>\mu_{j} = \frac{1}{l w_{j}} \sum_{i=1}^{l} g_{ij} x_{i}</math> | ||
+ | ##:: Вычислить <math>\Sigma_{j} = \frac{1}{l w_{j}} \sum_{i=1}^{l} g_{ij}(x_{i} - \mu_{j})(x_{i} - \mu_{j})^T</math> | ||
+ | ## Вычислить изменения логарифма правдоподобия <math>\Delta</math> | ||
+ | ## Если <math>\Delta < \varepsilon</math>, то досрочно выйти из цикла | ||
+ | # Вернуть <math>w_{j}, \mu_{j}, \Sigma_{j}, \ j = 1,\dotsc,k.</math> | ||
=== Последовательная сложность алгоритма === | === Последовательная сложность алгоритма === | ||
+ | Рассмотрим мультипликативную сложность одной итерации алгоритма: | ||
+ | * '''E-шаг:''' | ||
+ | ** Сложность обращения матрицы ковариаций и вычисления ее определителя — <math>O(n^{3})</math> (для простоты будем рассматривать метод Гаусса) | ||
+ | ** Сложность вычисления расстояния Махаланобиса (показателя экспоненты) при вычислении значения каждой скрытой переменной — <math>O(n^{2})</math> | ||
+ | ** Общая сложность E-шага — <math>O(k * n^{3} + k * l * n^{2})</math> | ||
+ | * '''M-шаг:''' | ||
+ | ** Сложность пересчета веса кластера — <math>O(l)</math> | ||
+ | ** Сложность пересчета центра кластера — <math>O(l * n)</math> | ||
+ | ** Сложность пересчета матрицы ковариаций кластера — <math>O(l * n)</math> | ||
+ | ** Общая сложность M-шага — <math>O(k * l * n)</math> | ||
+ | |||
+ | Таким образом, общая сложность одной итерации — <math>O(k * n^{2} * (l + n))</math>. При учете ограничения на максимальное число итераций получим оценку общей сложности алгоритма — <math>O(imax * k * n^{2} * (l + n))</math>. | ||
+ | |||
+ | Большое влияние на сложность E-шага оказывает необходимость обращать ковариационные матрицы и вычислять их определители. Помимо того, что это трудоемкая операция, ковариационные матрицы могут оказаться плохо обусловленными, что может привести к неустойчивости выборочных оценок параметров смеси. Обращения матриц можно избежать, если принять гипотезу о том, что в каждой компоненте смеси признаки некоррелированы, то есть ковариационные матрицы диагональные. В таком случае общая сложность алгоритма составит <math>O(imax * k * n * l)</math>. | ||
=== Информационный граф === | === Информационный граф === | ||
+ | На рисунках 2 и 3 представлены информационные графы<ref name=voevodin>Воеводин В.В., Воеводин Вл.В. Параллельные вычисления. - СПб.: БХВ-Петербург, 2002. - 608 с.</ref> шагов E и M соответственно. Прямоугольниками обозначены входные данные (с пометкой in/out), эллипсами обозначены операции, производимые над данными согласно алгоритму. | ||
+ | |||
+ | Как видно из рисунка 2, E-шаг распадается на: | ||
+ | * Обращение матриц ковариаций и вычисление их определителей | ||
+ | * Независимый расчет скрытых переменных по каждому объекту | ||
+ | |||
+ | [[File:E-step.png|thumb|center|1100px|Рисунок 2. Граф шага E с отображением входных и выходных данных.]] | ||
+ | |||
+ | На рисунке 3 видно, что некоторые промежуточные вычисления на шаге M можно проводить независимо для каждого объекта, но после этого придется проводить агрегацию полученных результатов по кластерам. | ||
+ | |||
+ | [[File:M-step.png|thumb|center|1100px|Рисунок 3. Граф шага M с отображением входных и выходных данных.]] | ||
=== Ресурс параллелизма алгоритма === | === Ресурс параллелизма алгоритма === | ||
+ | Из информационных графов на рисунках 2 и 3 видно, что: | ||
+ | * Шаг E — наиболее ресурсоемкий этап алгоритма — эффективно распараллеливается по объектам после предварительного обращения матриц ковариаций и вычисления их определителей: при расчете скрытых переменных для фиксированного объекта информация о других объектах не используется. Также стоит отметить, что для всех объектов выполняется одна и та же последовательность действий. В предположении о том, что число объектов выборки много больше размерности пространства признаков (что выполняется для подавляющего большинства прикладных задач), параллельная сложность E-шага — <math>O(k * n^{2})</math>. | ||
+ | * Менее ресурсоемкий шаг M также эффективно распараллеливается по объектам с последующей агрегацией результатов по кластерам: параллельная сложность M-шага — <math>O(k * n)</math>. | ||
+ | |||
+ | Таким образом, параллельная сложность EM-алгоритма — <math>O(imax * k * n^{2})</math>. | ||
=== Входные и выходные данные алгоритма === | === Входные и выходные данные алгоритма === | ||
+ | '''Входные данные:''' плотная матрица объектов-признаков <math>X \in \R^{l \times n}</math>, число кластеров <math>k \in \N</math>, максимальное число итераций <math>imax \in \N</math>, минимальная величина изменения логарифма правдоподобия <math>\varepsilon \in \R</math>. | ||
+ | |||
+ | '''Объём входных данных:''' <math>ln + 3</math> | ||
+ | |||
+ | '''Выходные данные:''' вещественный вектор весов <math>w \in \R^{k}</math>, <math>k</math> вещественных векторов математических ожиданий <math>\mu_{j} \in \R^{n}</math> и <math>k</math> вещественных ковариационных матриц <math>\Sigma_{j} \in \R^{n \times n}</math>. | ||
+ | |||
+ | '''Объём выходных данных:''' <math>k(n^{2}+n+1)</math> | ||
=== Свойства алгоритма === | === Свойства алгоритма === | ||
+ | Вычислительная мощность EM-алгоритма в предположении о том, что <math>l \gg n</math> (это условие выполняется в подавляющем большинстве прикладных задач), оценивается величиной <math>O(imax * k * n)</math>. | ||
+ | |||
+ | Достоинства EM-алгоритма: | ||
+ | * Мощная математическая основа | ||
+ | * Слабая чувствительность к выбросам | ||
+ | * Быстрая сходимость при удачной инициализации параметров | ||
+ | * Линейный рост сложности при увеличении количества объектов | ||
+ | |||
+ | Недостатки EM-алгоритма: | ||
+ | * В классическом варианте не может самостоятельно определить количество компонент смеси | ||
+ | * Трудоемкий и неустойчивый процесс обращения матриц ковариаций в случае несферических компонент смеси | ||
+ | * Не является детерминированным: в начале работы происходит инициализация начальных параметров случайным образом | ||
+ | * Не является устойчивым: результат работы сильно зависит от инициализации параметров. При неудачной инициализации может обладать низкой скоростью сходимости, может сойтись к локальному экстремуму или не сойтись вовсе | ||
== Программная реализация алгоритма == | == Программная реализация алгоритма == | ||
Строка 59: | Строка 156: | ||
=== Масштабируемость алгоритма и его реализации === | === Масштабируемость алгоритма и его реализации === | ||
+ | Исследование масштабируемости параллельной реализации EM-алгоритма проводилось на суперкомпьютере Blue Gene/P. Для этого использовалась собственная реализация алгоритма, написанная на языке C++ с использованием технологий MPI и OpenMP. | ||
+ | |||
+ | Для проведения исследования были сгенерированы наборы данных из смеси двух двумерных нормальных распределений с диагональными матрицами ковариаций. | ||
+ | |||
+ | Набор значений параметров запуска реализации алгоритма: | ||
+ | * число процессоров: 1, 4, 8, 16, 32, 48, 64, 80, 96, 112, 128 | ||
+ | * число элементов выборки: 1000, 25000, 50000, 75000, 100000, 200000, 300000, 400000, 500000, 600000, 700000, 800000, 900000, 1000000 | ||
+ | * размерность пространства признаков: 2 | ||
+ | * число кластеров: 2 | ||
+ | * максимальное число итераций: MAX_INT | ||
+ | * минимальная величина изменения логарифма правдоподобия: 1e-6 | ||
+ | |||
+ | Начальные значения параметров смеси задавались следующим образом: | ||
+ | * центры кластеров - случайные элементы выборки | ||
+ | * матрицы ковариаций - единичные | ||
+ | * веса компонент равны между собой | ||
+ | |||
+ | В связи с тем, что время и результат работы алгоритма сильно зависят от рандомизированной начальной инициализации центров кластеров, для каждого фиксированного набора параметров производилось по 25 запусков программы. В качестве итогового времени работы выбиралось минимальное время среди всех запусков. | ||
+ | |||
+ | На рисунке 4 изображена зависимость времени работы алгоритма от числа процессоров и числа элементов выборки. | ||
+ | |||
+ | [[File:EM-time.png|thumb|center|1200px|Рисунок 4. Зависимость времени работы параллельной реализации EM-алгоритма от числа процессоров и числа элементов выборки.]] | ||
+ | |||
+ | Как видно из рисунка 4, при фиксированном числе процессоров зависимость времени выполнения от количества объектов выборки очень близка к линейной, что хорошо согласуется с теоретической оценкой. | ||
+ | |||
+ | На следующих рисунках приведены графики зависимости производительности и эффективности реализации EM-алгоритма от числа процессоров и числа элементов выборки. | ||
+ | |||
+ | [[File:EM-performance.png|thumb|center|1200px|Рисунок 5. Зависимость производительности параллельной реализации EM-алгоритма от числа процессоров и числа элементов выборки.]] | ||
+ | |||
+ | [[File:EM-efficiency.png|thumb|center|1200px|Рисунок 6. Зависимость эффективности параллельной реализации EM-алгоритма от числа процессоров и числа элементов выборки.]] | ||
+ | |||
+ | В результате проведённых экспериментов был получен следующий диапазон эффективности реализации алгоритма: | ||
+ | * минимальная эффективность реализации 0.67% | ||
+ | * максимальная эффективность реализации 1.19% | ||
=== Динамические характеристики и эффективность реализации алгоритма === | === Динамические характеристики и эффективность реализации алгоритма === | ||
Строка 65: | Строка 196: | ||
=== Существующие реализации алгоритма === | === Существующие реализации алгоритма === | ||
+ | # [http://scikit-learn.org/stable/modules/generated/sklearn.mixture.GaussianMixture.html#sklearn.mixture.GaussianMixture scikit-learn(Python)] | ||
+ | # [http://docs.opencv.org/2.4/modules/ml/doc/expectation_maximization.html OpenCV(C++/Python)] | ||
+ | # [http://www.mlpack.org/docs/mlpack-git/doxygen.php?doc=classmlpack_1_1gmm_1_1GMM.html MLPack(C++)] | ||
+ | # [http://arma.sourceforge.net/docs.html#gmm_diag Armadillo(C++)] | ||
+ | # [https://www.mathworks.com/help/stats/fitgmdist.html Matlab] | ||
== Литература == | == Литература == |
Текущая версия на 15:49, 28 ноября 2016
Эта работа успешно выполнена Преподавателю: в основное пространство, в подстраницу Данное задание было проверено и зачтено. Проверено Dexter и Algoman. |
Авторы описания:
- Куваев А.С., группа 620 (математическая постановка задачи, схема реализации алгоритма, оценка сложности, описание ресурса параллелизма, реализация алгоритма)
- Щенявская Е.В., группа 616 (описание алгоритма, визуализация информационного графа и примера работы алгоритма, описание свойств алгоритма, замеры времени выполнения и визуализация результата)
Содержание
- 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 Общее описание алгоритма
Задача кластеризации заключается в разбиении входного множества объектов на непересекающиеся подмножества, называемые кластерами, так, чтобы каждый кластер состоял из схожих объектов, а объекты разных кластеров существенно отличались друг от друга[1].
Решение этой задачи принципиально неоднозначно по следующим причинам[1]:
- результат кластеризации зависит от способа задания меры сходства объектов выборки
- не существует однозначно наилучшего критерия качества кластеризации
- число кластеров, как правило, неизвестно заранее и задается из некоторых априорных соображений (хотя существуют алгоритмы, способные определять число кластеров автоматически)
По описанной выше причине существует большое число алгоритмов кластеризации, приводящих к различным разбиениям исходного множества объектов. На этой странице представлено описание одного из таких методов — EM-алгоритма. EM-алгоритм опирается на предположение о вероятностной природе данных: элементы выборки получены случайно и независимо из смеси распределений с фиксированным числом компонент [math]k[/math]. Таким образом, плотность распределения на множестве объектов имеет следующий вид:
- [math]p(x)=\sum_{j=1}^{k}w_{j}p_{j}(x), \ \sum_{j=1}^{k}w_{j}=1, \ w_{j} \ge 0 [/math], где [math]p_{j}[/math] - плотность распределения [math]j[/math]-й компоненты смеси (кластера).
Везде в дальнейшем будем предполагать, что [math]p_{j}[/math] имеют вид многомерных нормальных плотностей с произвольной матрицей ковариации: смеси нормальных распределений позволяют аппроксимировать произвольные непрерывные функции плотности с наперед заданной точностью[2]. Результатом работы EM-алгоритма являются оценки априорных вероятностей компонент смеси [math]w_{j}[/math], а также оценки векторов математических ожиданий и матриц ковариаций для каждой компоненты. Зная параметры распределения смеси, каждому объекту будет сопоставляться тот кластер, вероятность принадлежности к которому будет максимальной.
EM-алгоритм позволяет значительно упростить задачу максимизации правдоподобия выборки путем искусственного введения вспомогательной матрицы скрытых переменных [math]G[/math]. Алгоритм заключается в последовательном повторении шагов E(expectation) и M(maximization):
- На шаге E на основе текущего приближения параметров смеси по формуле Байеса вычисляются ожидаемые значения скрытых переменных [math]g_{ij}[/math] — апостериорные вероятности того, что [math]i[/math]-й объект принадлежит кластеру [math]j[/math].
- На шаге M решается задача максимизации правдоподобия для нахождения следующего приближения параметров смеси на основе текущего приближения и матрицы скрытых переменных, при этом решение этой задачи выписывается в явном виде.
1.2 Математическое описание алгоритма
Пусть заданы [math]l[/math] объектов [math]x_{1},\dotsc,x_{l}[/math], каждый из которых описывается [math]n[/math] числовыми признаками. Таким образом, определена матрица объектов-признаков [math]X \in \R^{l \times n}.[/math] Предполагается, что объекты выбраны случайно и независимо из смеси [math]n[/math]-мерных нормальных распределений с фиксированным числом компонент [math]k[/math]. Плотность распределения на множестве объектов имеет вид:
- [math]p(x)=\sum_{j=1}^{k}w_{j}p_{j}(x), \ \sum_{j=1}^{k}w_{j}=1, \ w_{j} \ge 0 [/math], где [math]p_{j}[/math] - плотность распределения [math]j[/math]-й компоненты смеси.
Каждая компонента смеси описывается [math]n[/math]-мерным вектором математических ожиданий [math]\mu_{j}[/math] и матрицей ковариаций [math]\Sigma_{j}[/math] порядка [math]n[/math], [math]j = 1,\dotsc,k[/math].
Входные данные: матрица объектов-признаков [math]X[/math], число кластеров [math]k[/math], максимальное число итераций [math]imax[/math], минимальная величина изменения логарифма правдоподобия [math]\varepsilon[/math].
Выходные данные: набор оценок весов, математических ожиданий и ковариационных матриц компонент смеси [math]\theta = (w_{1},...,w_{k}; \; \mu_{1},...,\mu_{k}; \; \Sigma_{1},...,\Sigma_{k})[/math], максимизирующий правдоподобие выборки.
Структура алгоритма:
- Инициализация параметров: существует большое число вариантов инициализации параметров распределения. Один из возможных подходов — задание векторов математических ожиданий компонент случайными элементами выборки, задание матриц ковариаций единичными матрицами и задание весов компонент равными [math]\frac{1}{k}.[/math]
- Последовательное выполнение шагов E и M до тех пор, пока правдоподобие выборки не стабилизируется или не будет достигнуто максимальное число итераций:
- Шаг E: вычисление значений скрытых переменных по формуле Байеса:
- [math]g_{ij} = \frac{w_{j} p_{j}(x_{i})}{\sum_{s=1}^{k} w_{s} p_{s}(x_{i})}[/math] , где [math]p_{j}(x_{i}) = \frac{1}{(2\pi)^{n/2}\sqrt{|\Sigma_{j}|}} \exp \biggl( -\frac{1}{2}(x_{i} - \mu_{j})^{T} \Sigma_{j}^{-1} (x_{i} - \mu_{j}) \biggr), \ i = 1,\dotsc,l; \ j = 1,\dotsc,k [/math]
- Шаг M: перерасчет параметров смеси на основе текущего приближения и матрицы скрытых переменных:
- [math]w_{j} = \frac{1}{l} \sum_{i=1}^{l} g_{ij}, \ j = 1,\dotsc,k;[/math]
- [math]\mu_{j} = \frac{1}{l w_{j}} \sum_{i=1}^{l} g_{ij} x_{i}, \ j = 1,\dotsc,k;[/math]
- [math]\Sigma_{j} = \frac{1}{l w_{j}} \sum_{i=1}^{l} g_{ij}(x_{i} - \mu_{j})(x_{i} - \mu_{j})^T, \ j = 1,\dotsc,k.[/math]
1.3 Вычислительное ядро алгоритма
Вычислительное ядро EM-алгоритма — процедура последовательного выполнения шагов E и M:
- На шаге E на основе текущего приближения параметров смеси вычисляются ожидаемые значения скрытых переменных
- На шаге M вычисляется следующее приближение параметров смеси на основе текущего приближения и матрицы скрытых переменных
Наиболее трудоемкой операцией с вычислительной точки зрения является шаг E, в ходе которого производится обращение ковариационных матриц, вычисление их определителей, а также многократное перемножение векторов и матриц при вычислении скрытых переменных.
1.4 Макроструктура алгоритма
Алгоритм состоит из двух итерационно проводимых макроопераций: E-шага и M-шага.
В ходе E-шага используется:
- Обращение ковариационной матрицы и вычисление ее определителя
- Умножение вектор-строки на матрицу и вектор-строки на вектор-столбец
В ходе M-шага используется:
- Умножение вектор-столбца на вектор-строку
- Взвешенное суммирование векторов и матриц
1.5 Схема реализации последовательного алгоритма
- Для всех [math]j = 1,\dotsc,k[/math]:
- Инициализировать [math]w_{j}, \ \mu_{j}, \ \Sigma_{j}[/math]
- Для всех [math]iter = 1,\dotsc,imax[/math]:
- Шаг E:
- Для всех [math]j = 1,\dotsc,k[/math]:
- Вычислить [math]|\Sigma_{j}|, \ \Sigma_{j}^{-1}[/math]
- Для всех [math]i = 1,\dotsc,l[/math]:
- [math]sum = 0[/math]
- Для всех [math]j = 1,\dotsc,k[/math]:
- Вычислить [math]p_{j}(x_{i}) = \frac{1}{(2\pi)^{n/2}\sqrt{|\Sigma_{j}|}} \exp \biggl( -\frac{1}{2}(x_{i} - \mu_{j})^{T} \Sigma_{j}^{-1} (x_{i} - \mu_{j}) \biggr)[/math]
- [math]sum = sum + w_{j}p_{j}(x_{i})[/math]
- Для всех [math]j = 1,\dotsc,k[/math]:
- Вычислить [math]g_{ij} = \frac{w_{j} p_{j}(x_{i})}{sum}[/math]
- Для всех [math]j = 1,\dotsc,k[/math]:
- Шаг M:
- Для всех [math]j = 1,\dotsc,k[/math]:
- Вычислить [math]w_{j} = \frac{1}{l} \sum_{i=1}^{l} g_{ij}[/math]
- Вычислить [math]\mu_{j} = \frac{1}{l w_{j}} \sum_{i=1}^{l} g_{ij} x_{i}[/math]
- Вычислить [math]\Sigma_{j} = \frac{1}{l w_{j}} \sum_{i=1}^{l} g_{ij}(x_{i} - \mu_{j})(x_{i} - \mu_{j})^T[/math]
- Для всех [math]j = 1,\dotsc,k[/math]:
- Вычислить изменения логарифма правдоподобия [math]\Delta[/math]
- Если [math]\Delta \lt \varepsilon[/math], то досрочно выйти из цикла
- Шаг E:
- Вернуть [math]w_{j}, \mu_{j}, \Sigma_{j}, \ j = 1,\dotsc,k.[/math]
1.6 Последовательная сложность алгоритма
Рассмотрим мультипликативную сложность одной итерации алгоритма:
- E-шаг:
- Сложность обращения матрицы ковариаций и вычисления ее определителя — [math]O(n^{3})[/math] (для простоты будем рассматривать метод Гаусса)
- Сложность вычисления расстояния Махаланобиса (показателя экспоненты) при вычислении значения каждой скрытой переменной — [math]O(n^{2})[/math]
- Общая сложность E-шага — [math]O(k * n^{3} + k * l * n^{2})[/math]
- M-шаг:
- Сложность пересчета веса кластера — [math]O(l)[/math]
- Сложность пересчета центра кластера — [math]O(l * n)[/math]
- Сложность пересчета матрицы ковариаций кластера — [math]O(l * n)[/math]
- Общая сложность M-шага — [math]O(k * l * n)[/math]
Таким образом, общая сложность одной итерации — [math]O(k * n^{2} * (l + n))[/math]. При учете ограничения на максимальное число итераций получим оценку общей сложности алгоритма — [math]O(imax * k * n^{2} * (l + n))[/math].
Большое влияние на сложность E-шага оказывает необходимость обращать ковариационные матрицы и вычислять их определители. Помимо того, что это трудоемкая операция, ковариационные матрицы могут оказаться плохо обусловленными, что может привести к неустойчивости выборочных оценок параметров смеси. Обращения матриц можно избежать, если принять гипотезу о том, что в каждой компоненте смеси признаки некоррелированы, то есть ковариационные матрицы диагональные. В таком случае общая сложность алгоритма составит [math]O(imax * k * n * l)[/math].
1.7 Информационный граф
На рисунках 2 и 3 представлены информационные графы[3] шагов E и M соответственно. Прямоугольниками обозначены входные данные (с пометкой in/out), эллипсами обозначены операции, производимые над данными согласно алгоритму.
Как видно из рисунка 2, E-шаг распадается на:
- Обращение матриц ковариаций и вычисление их определителей
- Независимый расчет скрытых переменных по каждому объекту
На рисунке 3 видно, что некоторые промежуточные вычисления на шаге M можно проводить независимо для каждого объекта, но после этого придется проводить агрегацию полученных результатов по кластерам.
1.8 Ресурс параллелизма алгоритма
Из информационных графов на рисунках 2 и 3 видно, что:
- Шаг E — наиболее ресурсоемкий этап алгоритма — эффективно распараллеливается по объектам после предварительного обращения матриц ковариаций и вычисления их определителей: при расчете скрытых переменных для фиксированного объекта информация о других объектах не используется. Также стоит отметить, что для всех объектов выполняется одна и та же последовательность действий. В предположении о том, что число объектов выборки много больше размерности пространства признаков (что выполняется для подавляющего большинства прикладных задач), параллельная сложность E-шага — [math]O(k * n^{2})[/math].
- Менее ресурсоемкий шаг M также эффективно распараллеливается по объектам с последующей агрегацией результатов по кластерам: параллельная сложность M-шага — [math]O(k * n)[/math].
Таким образом, параллельная сложность EM-алгоритма — [math]O(imax * k * n^{2})[/math].
1.9 Входные и выходные данные алгоритма
Входные данные: плотная матрица объектов-признаков [math]X \in \R^{l \times n}[/math], число кластеров [math]k \in \N[/math], максимальное число итераций [math]imax \in \N[/math], минимальная величина изменения логарифма правдоподобия [math]\varepsilon \in \R[/math].
Объём входных данных: [math]ln + 3[/math]
Выходные данные: вещественный вектор весов [math]w \in \R^{k}[/math], [math]k[/math] вещественных векторов математических ожиданий [math]\mu_{j} \in \R^{n}[/math] и [math]k[/math] вещественных ковариационных матриц [math]\Sigma_{j} \in \R^{n \times n}[/math].
Объём выходных данных: [math]k(n^{2}+n+1)[/math]
1.10 Свойства алгоритма
Вычислительная мощность EM-алгоритма в предположении о том, что [math]l \gg n[/math] (это условие выполняется в подавляющем большинстве прикладных задач), оценивается величиной [math]O(imax * k * n)[/math].
Достоинства EM-алгоритма:
- Мощная математическая основа
- Слабая чувствительность к выбросам
- Быстрая сходимость при удачной инициализации параметров
- Линейный рост сложности при увеличении количества объектов
Недостатки EM-алгоритма:
- В классическом варианте не может самостоятельно определить количество компонент смеси
- Трудоемкий и неустойчивый процесс обращения матриц ковариаций в случае несферических компонент смеси
- Не является детерминированным: в начале работы происходит инициализация начальных параметров случайным образом
- Не является устойчивым: результат работы сильно зависит от инициализации параметров. При неудачной инициализации может обладать низкой скоростью сходимости, может сойтись к локальному экстремуму или не сойтись вовсе
2 Программная реализация алгоритма
2.1 Особенности реализации последовательного алгоритма
2.2 Локальность данных и вычислений
2.3 Возможные способы и особенности параллельной реализации алгоритма
2.4 Масштабируемость алгоритма и его реализации
Исследование масштабируемости параллельной реализации EM-алгоритма проводилось на суперкомпьютере Blue Gene/P. Для этого использовалась собственная реализация алгоритма, написанная на языке C++ с использованием технологий MPI и OpenMP.
Для проведения исследования были сгенерированы наборы данных из смеси двух двумерных нормальных распределений с диагональными матрицами ковариаций.
Набор значений параметров запуска реализации алгоритма:
- число процессоров: 1, 4, 8, 16, 32, 48, 64, 80, 96, 112, 128
- число элементов выборки: 1000, 25000, 50000, 75000, 100000, 200000, 300000, 400000, 500000, 600000, 700000, 800000, 900000, 1000000
- размерность пространства признаков: 2
- число кластеров: 2
- максимальное число итераций: MAX_INT
- минимальная величина изменения логарифма правдоподобия: 1e-6
Начальные значения параметров смеси задавались следующим образом:
- центры кластеров - случайные элементы выборки
- матрицы ковариаций - единичные
- веса компонент равны между собой
В связи с тем, что время и результат работы алгоритма сильно зависят от рандомизированной начальной инициализации центров кластеров, для каждого фиксированного набора параметров производилось по 25 запусков программы. В качестве итогового времени работы выбиралось минимальное время среди всех запусков.
На рисунке 4 изображена зависимость времени работы алгоритма от числа процессоров и числа элементов выборки.
Как видно из рисунка 4, при фиксированном числе процессоров зависимость времени выполнения от количества объектов выборки очень близка к линейной, что хорошо согласуется с теоретической оценкой.
На следующих рисунках приведены графики зависимости производительности и эффективности реализации EM-алгоритма от числа процессоров и числа элементов выборки.
В результате проведённых экспериментов был получен следующий диапазон эффективности реализации алгоритма:
- минимальная эффективность реализации 0.67%
- максимальная эффективность реализации 1.19%
2.5 Динамические характеристики и эффективность реализации алгоритма
2.6 Выводы для классов архитектур
2.7 Существующие реализации алгоритма
3 Литература
- ↑ 1,0 1,1 Воронцов К.В., Математические методы обучения по прецедентам.
- ↑ Королёв В.Ю., ЕМ-алгоритм, его модификации и их применение к задаче разделения смесей вероятностных распределений. Теоретический обзор. - М.: ИПИРАН, 2007.
- ↑ Воеводин В.В., Воеводин Вл.В. Параллельные вычисления. - СПб.: БХВ-Петербург, 2002. - 608 с.