Эта статья посвящена кластеризации, а точнее, моему недавно добавленному в CRAN пакету ClusterR. Детали и примеры ниже в большинстве своем основаны на пакете Vignette.
Кластерный анализ или кластеризация — задача группирования набора объектов таким образом, чтобы объекты внутри одной группы (называемой кластером) были более похожи (в том или ином смысле) друг на друга, чем на объекты в других группах (кластерах). Это одна из главных задач исследовательского анализа данных и стандартная техника статистического анализа, применяемая в разных сферах, в т.ч. машинном обучении, распознавании образов, анализе изображений, поиске информации, биоинформатике, сжатии данных, компьютерной графике.
Наиболее известные примеры алгоритмов кластеризации — кластеризация на основе связности (иерархическая кластеризация), кластеризация на основе центров (метод 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)
Гауссова смесь распределений — статистическая модель для представления нормально распределенных субпопуляций внутри общей популяции. Гауссова смесь распределений параметризируется двумя типами значений — смесь весов компонентов и средних компонентов или ковариаций (для многомерного случая). Если количество компонентов известно, техника, чаще всего используемая для оценки параметров смеси распределений — ЕМ-алгоритм.
Функция GMM в пакете ClusterR — реализация на R класса для моделирования данных как гауссовой смеси распределений (GMM) из библиотеки Armadillo с предположением о диагональных ковариационных матрицах. Можно настраивать параметры функции, в том числе 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, здесь это справедливо для количества кластеров, равного 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-средних — метод векторного квантования, исходно применяемый в обработке сигналов, часто используется для кластерного анализа данных. Цель кластеризации методом k-средних — разделить n значений на k кластеров, в которых каждое значение принадлежит кластеру с ближайшим средним, выступающим прототипом кластера. Это приводит к разделению области данных на ячейки Вороного. Наиболее часто применяемый алгоритм использует итеративную уточняющую технику. Из-за повсеместного употребления, его называют алгоритмом k-средних; в частности, среди специалистов в сфере компьютерных наук он также известен как алгоритм Ллойда.
Пакет 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 на примере векторного квантования и пакета 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-средних в мини-группах — вариация классического алгоритма 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
- github.com/siddharth-agrawal/Mini-Batch-K-Means
- github.com/scikit-learn/scikit-learn/blob/51a765a/sklearn/cluster/k_means_.py#L1113
Автор: Инфопульс Украина
