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

EM-алгоритм
Последовательный алгоритм
Последовательная сложность	[math][/math]
Объём входных данных	[math]n * l + 1[/math]
Объём выходных данных	[math]n[/math]
Параллельный алгоритм
Высота ярусно-параллельной формы	[math][/math]
Ширина ярусно-параллельной формы	[math][/math]

Автор описания: Зинченко Д.А.

Алгоритм кластеризации, основанный на максимизации ожидания (EM-алгоритм)

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

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

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

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

1. Анализ неполных данных (данных с пропусками).

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

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

Дано множество объектов, каждый из которых представляет собой точку в n-мерном метрическом пространстве. Каждому измерению пространства соответствует некоторое свойство объекта. Необходимо разбить это множество на 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]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]\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 алгоритма с предварительной обработкой матриц ковариации:

// n - количество объектов

// m - количество координат объектов

// k - количество кластеров

Инициализировать веса и центры кластеров;
Инициализировать матрицы ковариации;

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

Наиболее часто используемым ограничением [3] является предположение, что матрица ковариации имеет диагональную форму. Тогда вычисление её определителя и обратной матрицы имеет сложнось [math]O(m)[/math], а суммарная сложность шага E и сложность всей итерации составит [math]O(k * m * n)[/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 Свойства алгоритма

EM-алгоритм обладает следующими преимуществами и недостатками:

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

1. Слабо чувствителен к выбросам

2. Прост в реализации

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

Недостатки:

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

2. Алгоритм не находит оптимальное количество кластеров. Количество кластеров, на которые нужно разбить множество объектов, является параметром алгоритма.

3. При большой размерности пространства объектов выдвинутое предположение о модели их распределения может быть некорректным.

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

2.1 Особенности реализации последовательного алгоритма

Последовательнй алгоритм может быть реализован следующим образом:

// n - количество объектов
// m - количество координат объектов
// k - количество кластеров

Инициализировать веса и центры кластеров

Инициализировать матрицы ковариации

// Основной цикл
while (условие остановки не выполнено) {
    // E-шаг
    for (l = 0; l < k; l++) {
        Вычислить определитель матрицы ковариации Sigma[l]
        Вычислить обратную матрицу матрице Sigma[l]  
    }
    
    for (i = 0; i < n; i++) {
        sum = 0;
             
        for (l = 0; l < k; l++) {
            Вычислить и запомнить p(x[i] | mu[l], Sigma[l]);
            sum += w[l] * p(x[i] | mu[l], Sigma[l])
        }
        
        for (l = 0; l < k; l++) {
            gamma[i, l] = w[l] * p(x[i] | mu[k], Sigma[k]) / sum;
        }
    }
    
    // M-шаг
    // обновление w[l]
    for (l = 0; i < k; l++) {
        cl_num[l] = 0;
        for (i = 0; i < n; i++) {
            cl_num[l] += gamma[i,l]
        }
        w[l] = cl_num[l] / n;
    }
    
    // обновление mu[l]
    for (l = 0; l < k; l++) {
        for (j = 0; j < m; j++) {
            mu[l][j] = 0;
            for (i = 0; i < n; i++) {
                 mu[l][j] += gamma[i,l] * x[i][j];
            }
            mu[l][j] /= cl_num[l];
        }
    }
    
    Обновить матрицы ковариации
}

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

Примеры реализаций последовательного алгоритма: OpenCV(C++) [4,5], Matlab [6], Scikit learn (Python) [7].

3 Литература

[1] 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] В.Ю.Королёв. ЕМ-алгоритм, его модификации и их применение к задаче разделения смесей вероятностных распределений.

[3] https://basegroup.ru/community/articles/em

[4] http://docs.opencv.org/2.4/modules/ml/doc/expectation_maximization.html

[5] https://github.com/opencv/opencv/blob/master/modules/ml/src/em.cpp

[6] http://www.mathworks.com/help/stats/clustering-using-gaussian-mixture-models.html

[7] http://scikit-learn.org/stable/modules/mixture.html

[8] Henggang Cui, Jinliang Wei, Wei Dai. Parallel Implementation of Expectation-Maximization for Fast Convergence - School of Computer Science, Carnegie Mellon University.