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

Файл:Symbol confirmed.svg

Эта работа успешно выполнена Преподавателю: в основное пространство, в подстраницу Данное задание было проверено и зачтено. Проверено Dexter и Algoman.

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

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

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

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

EM-алгоритм

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

где:

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

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

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

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

EM-алгоритм

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

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

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

     [math]\displaystyle{ \begin{align} & y_{ij} = \frac{w_j p_j(x_i)}{\sum_{l=1}^{k}w_lp_l(x_i)} \end{align} }[/math]

   • М-шаг: переоценка вектора параметров, используя текущие значения скрытых переменных:

     [math]\displaystyle{ \begin{align} & n_j = \sum_{i = 1}^{n}y_{ij}, \\ & \mu_j^{new} = \frac{1}{n_j} \sum_{i = 1}^{n}y_{ij}x_i, \\ & \Sigma_j^{new} = \frac{1}{n_j} \sum_{i = 1}^{n}y_{ij}(x_i - \mu_j^{new})(x_i - \mu_j^{new})^T,\\ & w_j^{new} = \frac{n_j}{n}. \end{align} }[/math]

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

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

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

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

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

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

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

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

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

Шаг E

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

Шаг M

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

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

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

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

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

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

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

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

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

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

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

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

E-шаг

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

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

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

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

Рисунок 2. Граф E-шага.

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

M-шаг

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

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

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

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

Рисунок 4. Граф M-шага.

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

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

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

Рисунок 7. Граф вычисления ковариационной матрицы [math]\displaystyle{ \Sigma_j }[/math].

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Литература

<references \>