EМ — масштабируемый алгоритм кластеризации

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

Его название происходит от слов «expectation-maximization», что можно перевести как «ожидание-максимизация». Это связано с тем, что каждая итерация содержит два шага – вычисление ожидаемых оценок параметров статистической модели (expectation) и максимизацию (maximisation) логарифмической функции правдоподобия для модели смеси статистических распределений, которая наилучшим образом аппроксимирует исходные данные. Алгоритм основан на методике итеративного вычисления оценок максимального правдоподобия, предложенной в 1977 г. (A. P. Demster, N. M. Laird, D. B. Rubin. Maximum Likelihood from Incomplete Data via the EM Algorithm).

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

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

Это значит, что EM-алгоритм даёт локально-оптимальное решение, что с точки зрения задачи кластеризации делает его эвристическим, т.е. дающим неточное и не единственное решение, но приемлемое в большинстве практически значимых случаев. Следовательно, формируемая с помощью алгоритма кластерная структура не является точной и единственно возможной, а является наиболее «правдоподобной». Поэтому подходить к её содержательной интерпретации следует с известной осторожностью.

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

Иными словами, предполагается, что данные в каждом кластере подчиняются определенному закону распределения, а именно, нормальному (рис. 1). С учетом этого предположения можно определить параметры распределения — математическое ожидание и дисперсию, которые соответствуют закону распределения значений признаков в кластере, наилучшим образом «подходящему» к исходным данным.

Таким образом, в основе алгоритма лежит предположение, что любое наблюдение исходного набора данных принадлежит ко всем формируемым кластерам, но с разной вероятностью. Тогда задача будет заключаться в «подгонке» распределений смеси к данным, а затем в определении вероятностей принадлежности наблюдения к каждому кластеру. Эти вероятности и будут результатом работы алгоритма. Очевидно, что наблюдение должно быть отнесено к тому кластеру, для которого данная вероятность выше.

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

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

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

Статистические основы алгоритма

Как отмечалось во введении, EM-алгоритм предполагает, что кластеризуемые данные подчиняются линейной комбинации (смеси) нормальных (гауссовых) распределений. Функция плотности вероятности нормального распределения непрерывной случайной величины x определяется по формуле:

f(x)=\frac{1}{^{\sqrt{2\cdot \pi \cdot \sigma ^{2}}}}\cdot exp\left \{ -\frac{(x-\mu_{x} )^{2}}{2\cdot \sigma _{x}^{2}} \right \}

где x∈(−∞;+∞), μ_xи \sigma _{x}^{2} — математическое ожидание и дисперсия случайной величины x соответственно.

Многомерное нормальное распределение для q-мерного пространства является обобщением предыдущего выражения. Многомерная функция плотности вероятности для случайного q-мерного вектора x=(x_1,x_2,...,x_q) имеет вид:

f(x)=\frac{1}{(2\cdot \pi)^{\frac{q}{2}}\left | \sum \right |^{\frac{1}{2}}}\cdot exp\left \{ -\frac{1}{2}\cdot (x-\mu )^{T}\cdot \sum ^{-1}\cdot (x-\mu )\right \}

где, μ — q-мерный вектор математических ожиданий x, \sum — его ковариационная матрица, \left | \sum \right | — определитель ковариационной матрицы, T — оператор транспонирования.

Введем в рассмотрение величину \delta ^{2}=(x-\mu )^{T}\cdot \left | \sum \right |^{-1}\cdot (x-\mu ), называемую квадратичным расстоянием Махаланобиса.

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

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

f(x)=\sum \limits_ {i=1}^{k}w_{i}\cdot f(x|i)

где f(x|i) — условное нормальное распределение для i-го кластера, w_i — доля наблюдений исходного набора данных, распределённых в i-й кластер (вес кластера).

Практическая реализация алгоритма

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

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

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

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

В процессе работы алгоритм вычисляет матрицу средних значений M, ковариаций R и весов W для функции распределения вероятности, описанной выше. Параметры, оцененные алгоритмом, сохраняются в таблице вида:

Следует отметить, что один из популярных алгоритмов кластеризации k-means можно рассматривать как частный случай EM-алгоритма, когда все элементы матрицы W одинаковы и w_{j}=\frac{1}{k}, а R=I, где I — единичная матрица (квадратная матрица, элементы главной диагонали которой равны единице, а остальные равны нулю).

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

Алгоритм начинает работу с инициализации, т.е. некоторого приближенного решения, которое может быть выбрано случайно или задано пользователем исходя из некоторых априорных сведений об исходных данных. Наиболее общим способом инициализации является присвоение элементам матрицы математических ожиданий M случайных значенийM\leftarrow \mu (т.е. случайно задаются центроиды кластеров), начальная ковариационная матрица определяется как единичная r\leftarrow I, веса кластеров задаются одинаковыми w_{j}=\frac{1}{k}.

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

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

INPUT:
   k // число кластеров
   j=1..k // номер кластера
   Y=\left \{ y_{1},y_{2},...,y_{n} \right \} // множество из n наблюдений q-мерного пространства
   i=1..n // номер наблюдения
   t=1..q// индекс измерения
   \varepsilon // минимальное приращение функции логарифмического правдоподобия по достижении которого алгоритм останавливается;
   Q_{max} // максимальное число итераций алгоритма (если алгоритм не сходится, то будет остановлен при достижении данного числа итераций).
OUTPUT:
   M, R, W // матрицы обновляемых параметров смеси распределений;
   X // матрица нормированных вероятностей принадлежности наблюдений к кластерам.
INITIALIZATION:
   M\leftarrow Random, R\leftarrow Random, W\leftarrow Random
WHILE
   \Delta llh\geq \varepsilon AND Q< Q_{max}
     PROCEDURE ЕXPECTATION
         llh=0
         FOR i=1 TO n
                sump=0
                FOR j=1 TO k
                       \delta _{ij}^{2}=(y_{i}-M_{j})^{T}\cdot R^{-1}\cdot (y_{i}-M_{j})
                      p_{ij}=\frac{w_{j}}{(2\cdot \pi)^{\frac{q}{2}}\cdot |R|^{\frac{1}{2}}}\cdot exp\left \{ -\frac{1}{2}\cdot \delta_{ij}^{2} \right \}
                      sump=sump+p_{ij}
                END FOR j
                X_{i}=\frac{P_i}{sump}
                \Delta llh=ln(sump)
                llh=llh+\Delta llh
                M=M+y_{i}\cdot X_{i}^{T}
                W=W+X_{i}
         END FOR i
  END PROCEDURE
 PROCEDURE MAXIMIZATION
         FOR j=1 TO k
               M_{j}=\frac{M_{j}}{W_{j}}
               FOR i=1 TO n
                       R_{j}=R_{j}+(y_{i}-M_{j})x_{ij}(y_{i}-M_{j})^{T}
               END FOR i
               R_{j}=\frac{R_{j}}{n}
               W=\frac{W}{n}
         END FOR j
  END PROCEDURE
END WHILE

Алгоритм содержит два шага: шаг ожидания (expectation) или Е-шаг и шаг максимизации (maximization) или M-шаг. Каждый из шагов повторяется до тех пор, пока изменение логарифмической функции правдоподобия \Delta llh не станет меньше, чем \varepsilon, или пока не будет достигнуто максимальное число итераций Q_{max}.

Принадлежность i-го наблюдения к j-му кластеру определяются условной вероятностью p_{ij}, причём количество этих вероятностей будет равно числу кластеров. Кластеры определяются своими центроидами и ковариационными матрицами. На основе центроида кластера и его ковариационной матрицы вычисляется расстояние Махаланобиса \delta_{ij}^{2} между вектором наблюдения y_i и центроидом кластера M_j. Для каждого кластера функция логарифмического правдоподобия llh определяется с помощью функции нормального распределения от расстояния Махаланобиса между центроидом и вектором наблюдения.

В процессе работы алгоритма итеративно модифицируются вероятности p_{ij}, центры кластеров M_j и ковариационные матрицы R_j. При этом на E-шаге перевычисляются вероятности, а на M-шаге - центроиды и ковариационные матрицы с использованием новых значений вероятности. После завершения работы алгоритма каждое наблюдение оказывается распределённым в кластер, для которого вероятность максимальна.

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

llh=\sum \limits_{i=1}^{n}ln\sum \limits_{j=1}^{k}p_{ij}

В псевдокоде сумма вероятностей по кластерам накапливается в переменной sump внутри цикла по j на E-шаге, а сумма по наблюдениям в цикле по i в переменной llh.

Буквами \delta ^{2}, R и Pобозначены матрицы квадратичных расстояний Махаланобиса, ковариаций и вероятностей принадлежности i-го наблюдения к j-му кластеру соответственно. Матрицы M, R и W являются вспомогательными временными матрицами и используются только для вычислений. При этом \left \| W \right \|=1, т.е. \sum \limits_{j=1}^{k} w_j (сумма весов всех кластеров в исходном наборе данных равна 1). Это означает, что любое случайно выбранное наблюдение принадлежит к какому-либо из кластеров с вероятностью 100%

Обозначение P_i, используемое в псевдокоде, представляет собой k-мерный вектор значений вероятностей p_{ij} принадлежности i-го наблюдения к j-му кластеру. Соответственно, x_{ij}обозначает нормированную вероятность принадлежности i-го наблюдения к j-му кластеру.

В матрице M размером q×k каждый элемент представляет собой оценку математического ожидания каждой входной переменной модели. А каждый столбец этой матрицы, вектор M_j — вектор оценок математического ожидания для всех переменных поj-му кластеру. R — диагональная матрица, у которой не равны нулю только элементы главной диагонали, т.е. R_{ij}=0 для всех i\neq j. Со статистической точки зрения это означает, что ковариации являются независимыми.

Именно диагональность ковариационной матрицы является ключевым предположением, которое делает алгоритм масштабируемым. В этом случае детерминант матрицы и его обращение может быть вычислено за время O(q), а алгоритм имеет сложность O(kqn). В случае недиагональной матрицы сложность алгоритма составит O(kq^{2}n), т.е. будет квадратично возрастать с увеличением размерности данных.

Важнейшей операцией, выполняемой на E-шаге алгоритма является вычисление квадратичных расстояний Махаланобиса \delta _{ij}^{2}. Если ковариационная матрица R является диагональной, то расстояние Махаланобиса от наблюдения y_i до центроида M_j(вектора средних по кластеру) j-го кластера с ковариационной матрицей R, будет:

\delta _{ij}^{2}=(y_{i}-M_{j})^{T}R^{-1}(y_{i}-M_{j})=\sum \limits_ {t=1}^{q}\frac{(y_{i}^{t}-M_{j}^{t})}{R_{tt}}

где t - индекс измерения q-мерного пространства, t=1..q, R^{-1}=\frac{1}{R}.

Если ковариационная матрица диагональная, то её обращение легко вычисляется, поскольку все её элементы вне главной диагонали равны 0. Ускорению вычислений также способствует то, что диагональная матрица может храниться в виде вектора её диагональных значений. Поскольку R не изменяется на Е-шаге, ее детерминант вычисляется только единожды, что делает вычисление вероятностей p_{ij} более быстрым. На M-шаге диагональность матрицы R также упрощает вычисления, поскольку недиагональные элементы матрицы (y_{i}-M_{j})x_{ij}(y_{i}-M_{j})^{T} равны нулю.

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

Численный эксперимент

Для иллюстрации работы EM-алгоритма и его сравнения с алгоритмом k-means рассмотрим результаты численного эксперимента. Для его проведения был использован набор данных, наблюдения которого представлены в виде точек в пространстве признаков на рис. 2.

Обратите внимание, что исходный набор данных не является простым для кластеризации, поскольку имеется явное перекрытие кластеров (области 1 и 2). В области 1 перекрываются кластеры 1 и 2, а в области 2 кластеры 4 и 5. Кластеры 3 и 6 расположены обособленно и, как ожидается, будут легко разделимыми. Для алгоритма k-means особые трудности должны возникнуть в местах перекрытия кластеров. Данное предположение подтверждается результатами, представленным на рис. 3.

В местах перекрытия кластеров наблюдается наибольшее число ошибок. В то же время обособленные кластеры 3 и 6 были распознаны алгоритмом k-means без ошибок. Как можно увидеть на рис. 4, алгоритм EM весьма успешно выявил перекрывающиеся кластеры, хотя и почти не распознал кластер 6.

Таким образом, можно сделать вывод, что алгоритм k-means может иметь преимущество при работе с обособленными (неперекрывающимися) кластерами, но полностью проигрывает алгоритму EM при наличии их перекрытия.

Литература

Другие материалы по теме:

Алгоритмы кластеризации на службе Data Mining

Кластеризация категорийных данных: масштабируемый алгоритм CLOPE