Участник:Александр Куваев/Алгоритм кластеризации, основанный на максимизации ожидания

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

Авторы описания:

• Куваев А.С., группа 620 (математическая постановка задачи, схема реализации алгоритма, оценка сложности, описание ресурса параллелизма, реализация алгоритма)
• Щенявская Е.В., группа 616 (описание алгоритма, визуализация информационного графа и примера работы алгоритма, описание свойств алгоритма, замеры времени выполнения и визуализация результата)

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

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

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

Решение этой задачи принципиально неоднозначно по следующим причинам:

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

По описанной выше причине существует большое число алгоритмов кластеризации, приводящих к различным разбиениям исходного множества объектов. На этой странице представлено описание одного из таких методов — 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] имеют вид многомерных нормальных плотностей с произвольной матрицей ковариации: смеси нормальных распределений позволяют аппроксимировать произвольные непрерывные функции плотности с наперед заданной точностью. Результатом работы 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. Визуализация компонент смеси после некоторых итераций алгоритма в двумерном пространстве признаков. Для каждой компоненты эллипсами выделены доверительные области с уровнями доверия 0.9, 0.95 и 0.99.

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 Схема реализации последовательного алгоритма

1. Для всех [math]j = 1,\dotsc,k[/math]:

   Инициализировать [math]w_{j}, \ \mu_{j}, \ \Sigma_{j}[/math]

2. Для всех [math]iter = 1,\dotsc,imax[/math]:
   1. Шаг 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]

   2. Шаг 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]

   3. Вычислить изменения логарифма правдоподобия [math]\Delta[/math]
   4. Если [math]\Delta \lt \varepsilon[/math], то досрочно выйти из цикла
3. Вернуть [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 представлены информационные графы шагов E и M соответственно. Прямоугольниками обозначены входные данные (с пометкой in/out), эллипсами обозначены операции, производимые над данными согласно алгоритму.

Как видно из рисунка 2, E-шаг распадается на:

• Обращение матриц ковариаций и вычисление их определителей
• Независимый расчет скрытых переменных по каждому объекту

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

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

Рисунок 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.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-алгоритма от числа процессоров и числа элементов выборки.

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

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

Рисунок 5. Зависимость производительности параллельной реализации EM-алгоритма от числа процессоров и числа элементов выборки.

Рисунок 6. Зависимость эффективности параллельной реализации EM-алгоритма от числа процессоров и числа элементов выборки.

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

• минимальная эффективность реализации 0.67%
• максимальная эффективность реализации 1.19%

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

1. scikit-learn(Python)
2. OpenCV(C++/Python)
3. MLPack(C++)
4. Armadillo(C++)
5. Matlab