Участник:Noite/EM-алгоритм кластеризации

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 Общее описание алгоритма

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

EM-алгоритм применяется для решения задач двух типов :

1. Анализ неполных данных (данных с пропусками).
2. Оценка максимального правдоподобия в случае, если функцию правдоподобия трудно исследовать аналитически. В этом случае введение набора скрытых переменных может существенно упростить задачу.

Задача кластеризации относится к задачам второго типа. В данной статье рассматривается следующая постановка задачи кластеризации:

Дано множество из N объектов, каждый из которых представляет собой точку в M-мерном метрическом пространстве. Каждому измерению пространства соответствует некоторое свойство объекта. Необходимо разбить это множество на K подмножеств (кластеров) так, чтобы, чтобы элементы каждого подмножества существенно отличались по свойствам от элементов других подмножеств.

Использование алгоритма EM подразумевает следующее предположение об объектах:

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

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

Итерация алгоритма состоит из двух последовательных шагов.

1. (Expectation) Вычисляются новые ожидаемые значения скрытых переменных.
2. (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]

Смесь распределений:

[math] P(x) = \sum_{j=1}^Kw_jp(x|\mu_{j}, \Sigma_{j})[/math]

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

[math]\ln P(x|w,\mu,\Sigma)= \sum_{n=1}^N \ln \{ \sum_{k=1}^K w_{k}p(x_{n}|\mu_{k},\Sigma_{k})\}[/math].

• В начале работы алгоритма задаются параметры начального приближения:

[math]w_{i} = 1/K, \mu_{i} = random(x_{j}), i=1..K, j=1..N[/math] [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]

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

• Далее итеративно выполняется следующая пара процедур:
  • 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]:

[math]\mu^{new}_{k}=\frac{1}{N_{k}}\sum_{n=1}^N \gamma_{n,k}x_{n}[/math], [math]\Sigma_{new}^{k}=\frac{1}{N_{k}}\sum_{n=1}^N \gamma_{n,k} (x_{n}-\mu_{k}^{new})(x_{n}-\mu_{k}^{new})^{T}[/math], [math]w_{k}^{new}=\frac{N_{k}}{N}[/math], [math]N_{k}=\sum_{n=1}^N \gamma_{n,k}[/math].

В EM алгоритме используется один из следующих критериев остановки:

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

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

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

• Обновление скрытых переменных на основе текущих параметров распределений
• Обновление параметров кластеров: весов, центров и матриц ковариации

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

Как записано в описании вычислительного ядра алгоритма, основную часть алгоритма составляют итеративно проводимые макрооперации: Expectation шаг - вычисление новых значений скрытых переменных, Maximization шаг - обновление параметров кластеров. Макроструктура алгоритма представлена на рисунке 1:

Рисунок 1: Макроструктура EM-алгоритма

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. Помимо этого сложность этих вычислений можно сократить еще сильнее, если ввести ограничения структуру матрицы ковариации [math]\Sigma[/math].

Наиболее часто используемым ограничением, позволяющим эффективно распараллеливать алгоритм, является предположение, что матрица ковариации имеет диагональную форму. Тогда вычисление её определителя и обратной матрицы имеет сложнось [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].

Рисунок 2: E-шаг EM алгоритма для одного объекта

На рисунке 3 представлена схема вычислений E-шага алгоритма для всех объектов: для объекта [math]x_{j}[/math] и параметров-векторов [math]w, \mu, \Sigma[/math] по формулам раздела 1.2 вычисляются значения скрытых переменных. [math]E_{i}, i=1..N[/math] - маркооперации, E-шаги для соответствующих объектов.

Рисунок 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].

Рисунок 4: M-шаг EM алгоритма для одного кластера

На рисунке 5 представлена схема вычислений M-шага алгоритма для всех кластеров. [math]M_{j}, j=1..K[/math] соответствует макрооперации - M-шагу для соответствующего кластера и скрытых переменных.

Рисунок 5: M-шаг EM алгоритма для всех кластеров

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-алгоритм обладает следующими преимуществами и недостатками:

Преимущества:

1. Слабо чувствителен к выбросам
2. Прост в реализации
3. Быстро сходится при удачном выборе начальных значений параметров.

Недостатки:

1. Неустойчив к выбору начальных значений параметров - от них зависит как скорость сходимости, так и результат работы: алгоритм может сойтись к локальному экстремуму функции правдоподобия, который может быть существенно ниже глобального.
2. Алгоритм не находит оптимальное количество кластеров. Количество кластеров, на которые нужно разбить множество объектов, является параметром алгоритма.
3. При большой размерности пространства объектов выдвинутое предположение о модели их распределения может быть некорректным.

1.10.2 Влияние структуры матрицы ковариации на форму кластера

Рисунок 6: Влияние структуры матрицы ковариации на форму кластеров для двумерного случая. a) - матрицы ковариации являются диагональными, b) - матрицы ковариации произвольны.

Рисунок 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 алгоритма:

1. центр распределения совпадает со случайно выбранным объектом [math]x_{i}[/math]
2. все наблюдения распределяются по кластерам случайным образом, а затем для каждого кластера пересчитывается цент соответствующего распределения
3. центры кластеров получаются в результате работы алгоритма k-means. Использование такой инициализации оправдано тем, что k-means требует меньших вычислительных затрат и быстрее сходится

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

Многие статьи, в которых представлены параллельные реализации EM алгоритма, описывают помимо использующихся программных и аппаратных решений результаты проведенных экспериментов с данными. В качестве примера в этом разделе будут рассмотрены результаты эксперементов, проведенных со связующим программным обеспечением (middleware) FREERIDE (FRamework for Rapid Implementation of Datamining Engines), обеспечивающим возможность распараллеливания как для систем с общей памятью (в частности, симметричных мультипроцессорных систем), так и для систем с распределенной памятью.

Эксперименты проводились на кластере машин Pentium 700MHz связанных при помощи Myrinet LANai 7.0 на трех наборах данных:

а) 3.3 * 10^6 объектов

b) 6.2 * 10^6 объектов

c) 12.5 * 10^6 объектов

Объекты представляют собой координаты точек 10-мерного пространства.

Количество машин в кластере варьировалось от 1 до 8, при этом на каждой машине запускалось до четырех параллельных потоков.

Результаты экспериментов представлены на рисунках 7-9:

Рисунок 7: набор данных a)

Рисунок 8: набор данных b)

Рисунок 9: набор данных c)

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

Данный раздел может быть дополнен

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

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

1. OpenCV(C++) (github)
2. Matlab
3. Scikit learn (Python)
4. MLPack(C++)

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

Здесь приводятся как конкретные реализации, так и статьи, описывающие подходы к реализациям распараллеливания EM-алгоритма

1. Apache Spark
2. Accord.NET Framework (C#)
3. Статья об использовании FREERIDE middleware для распараллеливания EM алгоритма
4. реализация с использованием CUDA
5. еще одна реализация с использованием CUDA
6. статья об использовании параллельного EM алгоритма в задачах компьютерного зрения
7. статья о распараллеливании EM алгоритма с использованием CUDA