Кластеризация с пакетом ClusterR, часть 1

Эта статья посвящена кластеризации, а точнее, моему недавно добавленному в CRAN пакету ClusterR. Детали и примеры ниже в большинстве своем основаны на пакете Vignette.

Кластерный анализ ^[1] или кластеризация — задача группирования набора объектов таким образом, чтобы объекты внутри одной группы (называемой кластером) были более похожи (в том или ином смысле) друг на друга, чем на объекты в других группах (кластерах). Это одна из главных задач исследовательского анализа данных и стандартная техника статистического анализа, применяемая в разных сферах, в т.ч. машинном обучении, распознавании образов, анализе изображений, поиске информации, биоинформатике, сжатии данных, компьютерной графике.

Наиболее известные примеры алгоритмов кластеризации — кластеризация на основе связности (иерархическая кластеризация), кластеризация на основе центров (метод k-средних, метод k-медоидов), кластеризация на основе распределений (GMM — Gaussian mixture models — Гауссова смесь распределений) и кластеризация на основе плотности (DBSCAN — Density-based spatial clustering of applications with noise — пространственная кластеризация приложений с шумом на основе плотности, OPTICS — Ordering points to identify the clustering structure — упорядочивание точек для определения структуры кластеризации, и др.).

Пакет ClusterR состоит из алгоритмов кластеризации на основе центров (метод k-средних, k-средних в мини-группах, k-медоидов) и распределений (GMM). Также пакет предлагает функции для:

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

Гауссова смесь распределений (GMM — Gaussian mixture models)

Гауссова смесь распределений ^[2] — статистическая модель для представления нормально распределенных субпопуляций внутри общей популяции. Гауссова смесь распределений параметризируется двумя типами значений — смесь весов компонентов и средних компонентов или ковариаций (для многомерного случая). Если количество компонентов известно, техника, чаще всего используемая для оценки параметров смеси распределений — ЕМ-алгоритм.

Функция GMM в пакете ClusterR — реализация на R класса для моделирования данных как гауссовой смеси распределений (GMM) из библиотеки Armadillo ^[3] с предположением о диагональных ковариационных матрицах. Можно настраивать параметры функции, в том числе gaussian_comps, dist_mode (eucl_dist, maha_dist), seed_mode (static_subset, random_subset, static_spread, random_spread), km_iter и em_iter (больше информации о параметрах в документации к пакету). Проиллюстрирую функцию GMM на синтетических данных dietary_survey_IBS.

library(ClusterR)
data(dietary_survey_IBS)
dim(dietary_survey_IBS)

## [1] 400  43

X = dietary_survey_IBS[, -ncol(dietary_survey_IBS)]   # данные (без зависимой переменной)
y = dietary_survey_IBS[, ncol(dietary_survey_IBS)]    # зависимая переменная
dat = center_scale(X, mean_center = T, sd_scale = T)  # центрирование и масштабирование данных
gmm = GMM(dat, 2, dist_mode = "maha_dist", seed_mode = "random_subset", km_iter = 10,
          em_iter = 10, verbose = F)          

# предсказание центроидов, ковариационной матрицы и весов
pr = predict_GMM(dat, gmm$centroids, gmm$covariance_matrices, gmm$weights)    

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

В дополнение к уже упомянутым функциям можно использовать Optimal_Clusters_GMM для оценки количества кластеров данных с помощью или информационного критерия Акаике (AIC, Akaike information), или байесовского информационного критерия (BIC, Bayesian information).

opt_gmm = Optimal_Clusters_GMM(dat, max_clusters = 10, criterion = "BIC", 
                               dist_mode = "maha_dist", seed_mode = "random_subset",
                               km_iter = 10, em_iter = 10, var_floor = 1e-10, 
                               plot_data = T)

При выборе модели из предопределенного множества предпочтительнее выбирать ту, у которой самый низкий BIC ^[4], здесь это справедливо для количества кластеров, равного 2.

Предполагая, что доступны метки истинности, для валидации выходных кластеров можно использовать методы external_validation (rand_index, adjusted_rand_index, jaccard_index, fowlkes_Mallows_index, mirkin_metric, purity, entropy, nmi (нормализованная взаимная информация) и var_info (вариация информации).

res = external_validation(dietary_survey_IBS$class, pr$cluster_labels, 
                          method = "adjusted_rand_index", summary_stats = T)
res

##  
## ---------------------------------------- 
## purity                         : 1 
## entropy                        : 0 
## normalized mutual information  : 1 
## variation of information       : 0 
## ---------------------------------------- 
## specificity                    : 1 
## sensitivity                    : 1 
## precision                      : 1 
## recall                         : 1 
## F-measure                      : 1 
## ---------------------------------------- 
## accuracy OR rand-index         : 1 
## adjusted-rand-index            : 1 
## jaccard-index                  : 1 
## fowlkes-mallows-index          : 1 
## mirkin-metric                  : 0 
## ----------------------------------------

И если параметр summary_stats установлен в TRUE, то возвращаются также метрики специфичности, чувствительности, точности, полноты, F-меры (specificity, sensitivity, precision, recall, F-measure соответственно).

Метод k-средних

Кластеризация методом k-средних ^[5] — метод векторного квантования, исходно применяемый в обработке сигналов, часто используется для кластерного анализа данных. Цель кластеризации методом k-средних — разделить n значений на k кластеров, в которых каждое значение принадлежит кластеру с ближайшим средним, выступающим прототипом кластера. Это приводит к разделению области данных на ячейки Вороного ^[6]. Наиболее часто применяемый алгоритм использует итеративную уточняющую технику. Из-за повсеместного употребления, его называют алгоритмом k-средних; в частности, среди специалистов в сфере компьютерных наук он также известен как алгоритм Ллойда ^[7].

Пакет ClusterR предоставляет две разных функции k-средних, KMeans_arma, реализацию на R метода k-средних из библиотеки armadillo, и KMeans_rcpp, которая использует пакет RcppArmadillo. Обе функции приводят к одним и тем же результатам, однако, возвращают разные признаки (код ниже иллюстрирует это).

KMeans_arma

KMeans_arma быстрее, чем функция KMeans_rcpp, однако, исходно она выводит центроиды только некоторых кластеров. Более того, количество колонок в данных должно превышать количество кластеров, иначе функция вернет ошибку. Кластеризация будет работать быстрее на многоядерных машинах с включенным OpenMP (например, -fopenmp в GCC). Алгоритм инициализируется один раз, и обычно 10 итераций достаточно для сходимости. Исходные центроиды распределяются с помощью одного из алгоритмов — keep_existing, static_subset, random_subset, static_spread или random_spread. Если seed_mode равно keep_existing, пользователь должен передать на вход матрицу центроидов.

Я уменьшу количество измерений в данных dietary_survey_IBS с помощью анализа главных компонент (PCA — principal component analysis), а именно — функции princomp из пакета stats, чтобы можно было вывести двумерный график кластеров, построенных в результате.

pca_dat = stats::princomp(dat)$scores[, 1:2]
km = KMeans_arma(pca_dat, clusters = 2, n_iter = 10, seed_mode = "random_subset", 
                 verbose = T, CENTROIDS = NULL)
pr = predict_KMeans(pca_dat, km)
table(dietary_survey_IBS$class, pr)
class(km) = 'matrix'
plot_2d(data = pca_dat, clusters = as.vector(pr), centroids_medoids = as.matrix(km)) 

KMeans_rcpp

Как утверждалось выше, функция KMeans_rcpp предлагает некоторые дополнительные возможности по сравнению с функцией KMeans_arma:

• она позволяет более одной инициализации (можно параллелизировать с OpenMP)
• помимо инициализаций с optimal_init, quantile_init, random или kmeans++, можно задать центроиды в параметре CENTROIDS
• время работы алгоритма и его сходимость можно настраивать параметрами num_init, max_iters и tol
• если num_init > 1, KMeans_rcpp возвращает атрибуты лучшей инициализации, используя критерий within-cluster-sum-of-squared-error (сумма-квадратов-ошибок-внутри-кластера)
• алгоритм возвращает следующие атрибуты: clusters, fuzzy_clusters (если fuzzy = TRUE), centroids, total_SSE, best_initialization, WCSS_per_cluster, obs_per_cluster, between.SS_DIV_total.SS

Больше подробностей о KMeans_rcpp есть в документации пакета. Проиллюстрирую разные параметры KMeans_rcpp на примере векторного квантования ^[8] и пакета OpenImageR.

library(OpenImageR)
path = 'elephant.jpg'
im = readImage(path)

# сначала изменим размер картинки, чтобы уменьшить количество измерений
im = resizeImage(im, 75, 75, method = 'bilinear')            
imageShow(im)                                                # вывод исходного изображения
im2 = apply(im, 3, as.vector)                                # векторизированный RGB 

# кластеризация с KMeans_rcpp
km_rc = KMeans_rcpp(im2, clusters = 5, num_init = 5, max_iters = 100, 
                    initializer = 'optimal_init', threads = 1, verbose = F)
km_rc$between.SS_DIV_total.SS

## [1] 0.9873009

Атрибут between.SS_DIV_total.SS равен (total_SSE — sum(WCSS_per_cluster)) / total_SSE. Если в кластеризации нет закономерности, промежуточная сумма квадратов будет очень небольшой частью общей суммы квадратов, а если атрибут between.SS_DIV_total.SS близок к 1.0, значения кластеризуются достаточно хорошо.

getcent = km_rc$centroids
getclust = km_rc$clusters
new_im = getcent[getclust, ]     # каждое значение ассоциируется с ближайшим центроидом
dim(new_im) = c(nrow(im), ncol(im), 3)        # обратное преобразование к трехмерной картинке
imageShow(new_im)

Кроме того, можно воспользоваться функцией Optimal_Clusters_KMeans (которая неявно использует KMeans_rcpp), чтобы определить оптимальное количество кластеров. Доступны такие критерии: variance_explained, WCSSE (within-cluster-sum-of-squared-error — сумма-квадратов-ошибок-внутри-кластера), dissimilarity, silhouette, distortion_fK, AIC, BIC и Adjusted_Rsquared. В документации по пакету есть больше информации по каждому критерию.

В следующем примере кода используется критерий distortion_fK, который полностью описан в статье “Selection of K in K-means clustering, Pham., Dimov., Nguyen., (2004)” («Выбор К в кластеризации методом К-средних»).

opt = Optimal_Clusters_KMeans(im2, max_clusters = 10, plot_clusters = T, 
                              verbose = F, criterion = 'distortion_fK', fK_threshold = 0.85)

Значения меньше фиксированной границы (здесь fK_threshold = 0.85) можно рекомендовать для кластеризации. Однако, существует более одного оптимального разбиения на кластеры, поэтому f(K) следует использовать только для предположения об их количестве, но окончательное решение, принимать это значение или нет — за пользователем.

Метод k-средних в мини-группах

Метод k-средних в мини-группах ^[9] — вариация классического алгоритма k-средних. Он особенно полезен для больших наборов данных, поскольку вместо использования всего набора (как в k-средних) берутся мини-группы из случайных значений данных, чтобы оптимизировать целевую функцию.

Параметры алгоритма MiniBatchKmeans почти такие же, как и функции KMeans_rcpp из пакета ClusterR. Наиболее важная разница — параметр batch_size (размер мини-групп) и init_fraction (процент данных для определения исходных центроидов, который применяется, если initializer равен ‘kmeans++’ или ‘quantile_init’).

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

path_d = 'dog.jpg'
im_d = readImage(path_d)

# сначала изменим размер картинки, чтобы уменьшить количество измерений
im_d = resizeImage(im_d, 350, 350, method = 'bilinear')            
imageShow(im_d)                                                # вывод исходного изображения

im3 = apply(im_d, 3, as.vector)                                # векторизированный RGB
dim(im3)                                              # исходные измерения данных

# 122500      3

Сначала выполним кластеризацию k-средних.

start = Sys.time()
km_init = KMeans_rcpp(im3, clusters = 5, num_init = 5, max_iters = 100, 
                    initializer = 'kmeans++', threads = 1, verbose = F)
end = Sys.time()
t = end - start
cat('time to complete :', t, attributes(t)$units, 'n')

# время выполнения : 2.44029 secs

getcent_init = km_init$centroids
getclust_init = km_init$clusters
new_im_init = getcent_init[getclust_init, ]  # каждое значение ассоциируется с ближайшим центроидом
dim(new_im_init) = c(nrow(im_d), ncol(im_d), 3)     # обратное преобразование к трехмерной картинке
imageShow(new_im_init)

Теперь кластеризацию k-средних в мини-группах.

start = Sys.time()
km_mb = MiniBatchKmeans(im3, clusters = 5, batch_size = 20, num_init = 5, max_iters = 100, 
                        init_fraction = 0.2, initializer = 'kmeans++', early_stop_iter = 10,
                        verbose = F)
pr_mb = predict_MBatchKMeans(im3, km_mb$centroids)
end = Sys.time()
t = end - start
cat('time to complete :', t, attributes(t)$units, 'n')

# время выполнения : 0.8346727 secs 

getcent_mb = km_mb$centroids
new_im_mb = getcent_mb[pr_mb, ]   # каждое значение ассоциируется с ближайшим центроидом
dim(new_im_mb) = c(nrow(im_d), ncol(im_d), 3)     # обратное преобразование к трехмерной картинке
imageShow(new_im_mb)

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

Чтобы реализовать метод k-средних в мини-группах, я использовал следующие ресурсы:

• www.eecs.tufts.edu/~dsculley/papers/fastkmeans.pdf ^[9]
• github.com/siddharth-agrawal/Mini-Batch-K-Means ^[10]
• github.com/scikit-learn/scikit-learn/blob/51a765a/sklearn/cluster/k_means_.py#L1113 ^[11]

Автор: Инфопульс Украина

Источник ^[12]