CC BY
46
№ 4 (52) 2014
М. А. Баранов, аспирант Национального исследовательского университета «Высшая школа экономики», г. Москва, thenorthcat@gmail.com
Параллельная версия жадного алгоритма кластеризации
В статье предлагается параллельная версия одного из алгоритмов кластеризации, принцип работы которого основан на так называемом жадном подходе . Для решения задачи распараллеливания алгоритма была выбрана технология CUDA, разработанная компанией NVIDIA . Приводятся программный код и результаты вычислительных экспериментов для матриц схожести разного размера
Ключевые слова: кластеризация, жадный алгоритм, параллельные вычисления .
введение
Одной из наиболее важных задач в информационном поиске является кластеризация — разбиение исходного множества объектов на группы, состоящие из схожих объектов. Кластеризация нашла широкое применение в различных областях знаний: в биологии, социологии, информатике, астрономии, медицине, археологии, маркетинговых исследованиях.
Кластерный анализ относится к процессам обучения без учителя, его основная цель — разбить множество объектов на группы таким образом, чтобы объекты внутри одной группы были максимально похожими друг на друга, но в то же время максимально отличались от объектов другой группы.
В настоящее время разработано множество алгоритмов кластеризации, использующих различные подходы к решению задачи кластерного анализа. Их классификация подробно изложена в работе [4]. Из всего многообразия используемых при кластеризации подходов стоит выделить так называемые жадные алгоритмы, суть работы которых сводится к тому, что на каждом шаге они делают локально оптимальный выбор в расчете на то, что это приведет к оптимальному решению всей задачи [5]. Жад-
ные алгоритмы часто используются при решении задач кластеризации [7, 10].
Одним из таких алгоритмов является алгоритм, предложенный в работах [1, 2] (далее он будет именоваться Greedy). Там же показана эффективность данного алгоритма при кластеризации коллекций текстовых документов. На вход алгоритма подаются матрица схожести документов и пороговое значение степени схожести (параметр threshold). Блок-схема алгоритма представлена на рис. 1.
Целью данной работы является сравнительный анализ параметров параллельной и последовательной версий алгоритма Greedy при кластеризации коллекций документов разумного размера (несколько тысяч документов).
краткое описание технологии CUDA
Для решения задачи распараллеливания алгоритма была выбрана аппаратно-программная архитектура CUDA (Compute Unified Device Architecture), разработанная компанией NVIDIA. Отметим, ранее в ряде работ данная технология была успешно использована для создания параллельных версий алгоритмов кластеризации [8, 9].
Архитектура CUDA основана на архитектуре вычислительных систем SIMD (single
№ 4 (52) 2014
Да
Добавить документ в кластер
Нет
> гДа
Вернуть кластерное разбиение
1
КОНЕЦ
Рис. 1. Блок-схема алгоритма Greedy
59
-n journal of applied informatics
№ 4 (52) 2014 ' -
1 2 3 4
5 6 7 8
9 10 11 12
=>
1 2 3 4 5 6 7 8 9 10 11 12
Рис. 2. Упаковка матрицы в одномерный массив
instruction, multiple data — одиночный поток команд и множественный поток данных), по классификации Флинна.
Определим основные термины, применяемые в технологии CUDA. Центральный процессор, управляющий выполнением программы, называется хостом (host). Видеокарта, выступающая в роли сопроцессора, называется устройством (device). Набор операций, выполняемых в определенном порядке, называется потоком (thread). Потоки объединены в блоки (block), которые, в свою очередь, объединены в грид (grid). Функция, выполняемая на устройстве, называется ядром (kernel).
Отличие потоков графического процессора GPU (Graphics Processing Unit) от потоков CPU состоит в следующем: 1) потоки GPU характеризуются очень низкими накладными расходами на создание; 2) переключение контекста потоков GPU происходит гораздо быстрее, чем потоков CPU; 3) GPU исполняет множество потоков одновременно.
Параллельные части алгоритма выполняются непосредственно на устройстве. Центральный процессор отвечает за выполнение последовательных операций, таких как выделение памяти, копирование данных, запуск ядер.
Таким образом, типичная программа на CUDA состоит из следующих операций:
• выделение памяти на устройстве с помощью функции cudaMalloc;
• копирование данных на устройство с помощью функции cudaMemcpy;
• запуск ядер для обработки данных;
• копирование результатов вычислений с устройства с помощью функции cudaMemcpy.
Приложения для CUDA разрабатываются на упрощенном диалекте языка C. Для трансляции используется компилятор nvcc, входящий в состав CUDA SDK.
Параллельная версия алгоритма Greedy
Код главного алгоритма Greedy представлен в листинге 1.
На вход функции greedy_clustering подаются два одномерных массива matrix и clusters и переменная MATRIX_SIZE. Массив matrix содержит строки матрицы схожести (рис. 2), массив clusters используется для сохранения результатов кластеризации, переменная MATRIX_SIZE задает число элементов массива, в который отображается матрица схожести. Число элементов в этом массиве соответствует числу документов в текстовой коллекции, i-й индекс массива соответствует i-му документу коллекции, а значение элемента массива с этим индексом является номером кластера, которому принадлежит документ. Изначально данный массив должен быть заполнен нулями.
Данная параллельная реализация требует, чтобы значение меры схожести лежало в диапазоне от 0 до 100, где 0 означает отсутствие какой-либо схожести документов, а 100 — полное сходство.
Дадим краткие пояснения к листингу 1. В строках 3 и 4 устанавливается количество потоков в блоке и вычисляется необходимое число блоков. Поскольку максимальное количество потоков в блоке определяется конкретным устройством1, данный параметр выбирался таким образом, чтобы программный
1 http://en.wikipedia.org/wiki/CUDA#Version_features. and_specifications.
№ 4 (52) 2014
Листинг 1
1. void greedy_clustering( unsigned char *matrix, unsigned int *clusters ){
2.
3. const int THREADS = 512;
4. const int BLOCKS = MATRIX_SIZE/THREADS + 1;
5.
6. unsigned char *d_matrix;
7. unsigned int*d_clusters;
8. unsigned int*d_cluster_id;
9. unsigned int*d_sum_vec;
10. unsigned int*d_stencil;
11.
12. unsigned int cluster_id = 1;
13.
14. cudaMalloc((void**)&d_matrix, MATRIX_SIZE*
15. MATRIX_SIZE*sizeof(unsigned char) ) ;
16. cudaMalloc((void**)&d_clusters,
17. MATRIX_SIZE*sizeof(unsigned int) );
18. cudaMalloc((void**)&d_cluster_id, sizeof(unsigned int) ) ;
19. cudaMalloc((void**)&d_sum_vec, MATRIX_SIZE*
20. sizeof(unsigned int) );
21. cudaMalloc((void**)&d_stencil, MATRIX_SIZE*sizeof(unsigned int));
22.
23. cudaMemcpy( d_matrix, matrix, MATRIX_SIZE*MATRIX_SIZE
24. *sizeof(unsigned char), cudaMemcpyHostToDevice ) ;
25. cudaMemcpy( d_clusters, clusters, MATRIX_SIZE*sizeof(unsigned int),
26. cudaMemcpyHostToDevice ) ;
27.
28. unsigned int min_value = 0;
29.
30. find_rows_sum<<<BLOCKS,THREADS>>>( MATRIX_SIZE,
31. d_matrix, d_sum_vec );
32. set_stencil<<<BLOCKS,THREADS>>>>(d_stencil);
33.
34. thrust::device_vector<unsigned int> filtered_vector(MATRIX_SIZE);
35. thrust::fill(filtered_vector.begin(), filtered_vector.end(), (int) 0);
36.
37. while ( min_value == 0 ){
38. unsigned int index = max_index( d_sum_vec );
39.
40. greedy<<<BLOCKS,THREADS>>>( d_matrix, index, d_clusters,
41. cluster_id, THRESHOLD, MATRIX_SIZE );
42.
43. get_doc_indexes_in_current_cluster( d_clusters ,cluster_id, d_stencil, filtered_vector);
44. unsigned int* indices_array =
45. thrust::raw_pointer_cast( &filtered_vector[0] );
46. reduce_sum_vec(indices_array, d_matrix, d_sum_vec);
47.
48. zero_matrix<<<BLOCKS,THREADS>>>( d_matrix, d_clusters,
49. MATRIX_SIZE );
50.
51. min_value =find_min_value( d_clusters, MATRIX_SIZE );
61
№ 4 (52) 2014
Окончание листинга 1
52. ++cluster_ id;
53.
54. thrust::fill(filtered_vector.begin(), filtered_vector.end(), (int) 0);
55. }
56.
57.
58. cudaMemcpy( clusters, d_clusters, MATRIX_SIZE*sizeof(unsigned int),
59. cudaMemcpyDeviceToHost ) ;
60.
61. cudaFree( d_matrix ) ;
62. cudaFree( d_clusters );
63. cudaFree( d_cluster_id ) ;
64. cudaFree( d_sum_vec );
65. cudaFree( d_stencil );
66. }
код мог быть исполнен на любой видеокарте, поддерживающей технологию CUDA.
В строках 14-21 выделяется память под переменные и массивы, которые будут использоваться на устройстве. В строках 23 - 26 входные данные копируются в память устройства.
Цикл while (строки 37-55) повторяется до тех пор, пока в матрице не останется элементов, не включенных в какой-либо кластер (т. е. в массиве d_clusters не останется нулевых элементов).
* * * *
Функция find_rows_sum (рис. 3 и листинг 2) суммирует значения элементов в каждой строке матрицы схожести и сохраняет полученный результат в массиве d_sum_vec. При этом каждый поток обрабатывает свою строку матрицы схожести.
Функция set_stencil (рис. 4 и листинг 3) в каждый элемент массива заносит его индекс.
Полученный результат используется затем в функции get_doc_indexes_in_current_ cluster.
Поток 1 * Поток 2 * Поток 3 * Поток 4 *
1 0,6 0,7 0,8
0,6 1 0,9 0,3
0,7 0,9 1 0,5
0,8 0,3 0,5 1
3,1
2,8
3,1
2,6
d_matrix find_rows_sum() d_sum_vce
Рис. 3. Работа функции find_rows_sum
62
^^ I 0 I 1 2 I 3 4 5 I 6 7 8
01 2345678
d_stencil.
Внизу - индексы массива set_stencil()
d_stencil
Рис. 4. Работа функции set_stencil
№ 4 (52) 2014
Функция maxjndex (листинг 4) находит индекс максимального элемента в массиве d_sum_vec, т. е. фактически индекс строки в матрице схожести с максимальной суммой элементов.
Входом функции greedy (листинг 5) служат матрица (d_matrix), индекс максимальной строки (index), массив с результатами кластеризации (d_clusters) и номер текущего кластера (clusterjd). Каждый поток обрабатывает свой документ в строке. Если
степень схожести документа с центром кластера больше значения параметра threshold, то соответствующему этому документу элементу массива d_clusters присваивается значение параметра duster_id.
Функция get_doc_indexes_in_current_ cluster (рис. 5) находит индексы документов, включенных в кластер с заданным номером. Результат заносится в динамический массив device_vector. В массиве d_stencil содержатся индексы документов. Исполь-
Листинг2
1 2
3
4
5
6
7
8
9
10 11 12
13
14
15
16
17
18
_global_ void find_rows_sum( unsigned int MATRIX_SIZE,
unsigned char *d_matrix, unsigned int *d_sum_vec){
int threadNumber = threadIdx.x + blockDim.x*blockIdx.x;
if ( threadNumber >= MATRIX_SIZE ){
return;
}
long int sum = 0;
for ( int column = 0; column < MATRIX_SIZE; ++column ){ sum += d_matrix[ column + threadNumber*MATRIX_SIZE ];
}
d_sum_vec[ threadNumber ] = sum; return;
_global_ void set_stencil(int*d_stencil){
unsigned int threadNumber = threadIdx.x + blockDim.x*blockIdx.x; d_stencil[threadNumber] = threadNumber;
}
1. unsigned int max_index( unsigned int *d_sum_vec){
2.
thrust::device_ptr<unsigned int> d_ptr(d_sum_vec);
3
4
5
6
7
8 9
10. unsigned int position = iter - iVec.begin();
11.
12. return position;
13. }
thrust::device_vector<unsigned int> iVec(d_ptr, d_ptr +
MATRIX_SIZE); thrust::device_vector<int>::iterator iter =
thrust::max_element(iVec.begin(), iVec.end());
Листинг3
Листинг 4
№ 4 (52) 2014
Листинг 5
1. _global_ void greedy (unsigned char *d_matrix, int index,
2. unsigned int *d_clusters, unsigned int cluster_id,
3. unsigned char THRESHOLD, unsigned int MATRIX_SIZE){
4.
5. int threadNumber = threadIdx.x + blockDim.x*blockIdx.x;
6.
7. if ( threadNumber >= MATRIX_SIZE ){
8. return;
9. }
10.
11. if ( d_matrix[ threadNumber + index*MATRIX_SIZE ]
12. > THRESHOLD )
13. {
14. d_clusters[ threadNumber ] = cluster_id;
15. }
16. }
0 1 2 1 3 1 2 0 1
d clusters
0 1 2 3 4 5 6 7 8
d stencil
1 3 5 8
пусто
>
get_doc_indexes_in_current_cluster()
filtered vector
filtered vector
<
cluster_id =1
Рис. 5. Работа функции get_doc_indexes_in_current_cluster
зование этого массива необходимо, так как функция copy_if из библиотеки thrust копирует значения элементов, а не их индексы.
Программный код функции get_doc_ indexes_in_current_cluster приведен на листинге 6.
В качестве последнего параметра в функцию copy_if передается предикат, управляющий процессом фильтрации, в виде следующей структуры (листинг 7).
Таким образом, элемент из массива d_stencil включается в массив filtered_vector при условии равенства этого элемента значению clusterjd.
Элементы, соответствующие включенным в кластер документам, должны быть обнулены в матрице схожести. Затем необходимо вновь просуммировать значения элементов в каждой строке матрицы и сохранить по-
лученный результат в массиве d_sum_vec. Чтобы избежать лишних операций, функция reduce_sum_vec (листинг 8 и рис. 6) уменьшает соответствующие суммы для каждой строки непосредственно в массиве d_sum_ vec без повторного вычисления сумм элементов.
В функции zero_matrix (см. листинг 9 и рис. 6) в матрице схожести (параметр d_matrix) обнуляются элементы, соответствующие включенным в кластер документам. При этом каждый поток обрабатывает только один определенный документ (номер обрабатываемого документа соответствует номеру потока threadNumber).
Функция find_min_value (см. листинг 10) находит минимальное значение в массиве d_clusters, содержащем результаты кластеризации. Первоначально данный массив за-
№ 4 (52) 2014
Листинг 6
void get_doc_indexes_in_current_cluster( unsigned int * d_clusters, unsigned int cluster_id, unsigned int *d_stencil, thrust::device_vector<unsigned int> &filtered_vector)
{
thrust::device_ptr<int> dev_ptr = thrust::device_pointer_cast(d_clusters); thrust::device_ptr<int> stencil_ptr =
thrust::device_pointer_cast(d_stencil);
my_equal e;
e.val = cluster_id;
thrust::copy_if(stencil_ptr, stencil_ptr+MATRIX_SIZE, dev_ptr, filtered_vector.begin(), e);
}
Листинг 7
struct my_equal{
int val;
bool operator() ( const int x)
{
return (x = = val);
}
};
Листинг 8
1. _global_ void reduce_sum_vec(unsigned int *d_doc_indexes, unsigned
2. char *d_matrix, unsigned int * d_sum_vec, unsigned int MATRIX_SIZE)
3. {
4. int threadNumber = threadIdx.x + blockDim.x*blockIdx.x;
5.
6. if ( threadNumber >= MATRIX_SIZE ){
7. return;
8. }
9.
10. for (int i = 0; i<MATRIX_SIZE; ++i){
11. if ( (d_doc_indexes[i] == 0) && (i>0) ){
12. return;
13. }
14. int el = d_matrix[d_doc_indexes[i] + threadNumber*MATRIX_SIZE];
15. if (d_sum_vec[ threadNumber ] > 0){
16. d_sum_vec[ threadNumber ] -= el;
17. }
18. if(d_doc_indexes[i] == threadNumber){
19. d_sum_vec[ threadNumber ] = 0;
20. } 21. }
22. }
№ 4 (52) 2014
d_matnx d_sum_vec
Рис. 6. Работа функций zero_matrix и reduce_sum_vec
Листинг 9
1 0,6 0,7 0,8
0,6 1 0,8 0,3
0,7 0,8 1 0,5
0,8 0,3 0,5 1
d matrix
3,1
2, 7
3, 0
2,6
d sum vec
В кластер включены документы 1 и 4 zero_matrix() reduce_sum_vec()
0 0 0 0
0 1 0,8 0
0 0,8 1 0
0 0 0 0
0
1,8
1,8
0
полнен нулями. Если в массиве d_clusters не остается ни одного нулевого элемента, это означает, что все документы были распределены по кластерам, и процесс кластеризации следует завершить.
вычислительные эксперименты
Вычислительные эксперименты проводились с использованием видеокарты Tesla C2050 (табл. 1) и процессора Intel Core i7-
№ 4 (52) 2014
Листинг 10
unsigned int find_min_value( unsigned int *d_clusters, unsigned int MATRIX_SIZE ){
thrust::device_ptr<unsigned int> d_ptr(d_clusters);
thrust::device_vector<unsigned int> iVec(d_ptr, d_ptr + MATRIX_SIZE);
thrust::device_vector<unsigned int>::iterator iter =
thrust::min_element(iVec.begin(), iVec.end());
unsigned int min_value return min_value;
iter;
Таблица 1
Таблица 2
Основные характеристики GPU NVIDIA Tesla C20502 3
Характеристика Значение
Число ядер CUDA 448
Частота работы ядер CUDA, ГГц 1,15
Полный объем специальной памяти, Гб 3
Тактовая частота графического процессора, ГГц 1,5
Количество GPU 14
Максимальное число потоков в блоке 1024
Максимальная размерность блока потоков (x, У, z) 1024 х 1024 х 64
Максимальная размерность грида из блоков 65535 х 65535 х 65535
Год выпуска 2010
3517и (табл. 2). Часть кода, связанная с обработкой ошибок, в листингах программ с целью уменьшения их объема опущена.
Основные характеристики CPU Intel Core i7-3517U4
Характеристика Значение
Количество ядер 2
Количество потоков 4
Тактовая частота процессора, МГц 1900-3000
Кэш 1-го уровня, Кбайт 128
Кэш 2-го уровня, Кбайт 512
Кэш 3-го уровня, Кбайт 4096
Год выпуска 2012
В таблице 3 приводятся результаты измерения времени выполнения последовательной (реализованной на С++) и параллельной (CUDA) версий алгоритма Greedy. Для параллельной версии приведены данные по времени выполнения как с учетом операций копирования данных с хоста на устройство и обратно, так и без учета этих операций.
Для наглядности на рис. 7 приведены графики зависимости времени выполнения параллельной и последовательной версий алгоритма Greedy от размера коллекции.
2 http://en.wikipedia.org/wiki/CUDA.
3 http://www.nvidia.ru/object/product_tesla_C2050_
C2070 ru.html.
4 http://ark.intel.com/ru/products/65714/Intel-Core-i7-3517U-Processor-4M-Cache-up-to-3_00-GHz.
=v 67
Таблица 3
Оценка времени выполнения алгоритма Greedy (реализация на C++ и CUDA)
Размер матрицы схожести, N х N Среднее время выполнения, мс
С++ CUDA, с учетом операций копирования CUDA, без учета операций копирования
100х 100 < 1 86 3
500х500 < 1 155 7
1000х 1000 4 206 11
5000х5000 386 286 69
10000х 10000 1 738 491 339
15000х 15000 4 045 1 054 943
20000х20000 7 769 1 495 1 379
25000х25000 15 214 2 926 2 453
30 000 х 30 000 23 439 4 093 3 678
25000
C++
CUDA (с копированием) CUDA (без копирования)
5000
10000
15000
20000
25000
30000
Число документов, шт.
Рис. 7. График зависимости времени выполнения параллельной и последовательной версий
алгоритма Greedy от размера коллекции
Заключение
При относительно большом размере матрицы схожести параллельная версия алгоритма Greedy работает быстрее последовательной версии, и различие во вре-
мени выполнения тем больше, чем больше размер матрицы схожести. Однако при небольших размерах матрицы схожести более быстрой оказывается последовательная версия. Такой результат можно объяс-
0
№ 4 (52) 2014
нить дополнительными расходами времени на копирование данных на устройство CUDA и обратно, а также тем, что отдельно взятое ядро графического процессора значительно уступает ядру CPU в скорости работы.
Таким образом, при кластеризации коллекций, состоящих из более чем 5000 документов, целесообразно использовать параллельную версию алгоритма Greedy.
Отметим, что параллельная версия работает примерно в 6 раз быстрее последовательной. Поскольку в экспериментах использовались процессор 2012 г. выпуска и видеокарта 2010 г. выпуска, то на более современной видеокарте следует ожидать увеличения скорости выполнения.
Список литературы
1. Барахнин В. Б., Нехаева В. А., Федотов А. М. О задании меры сходства для кластеризации текстовых документов // Вестник Новосибирского государственного университета. Серия: Информационные технологии. 2008. Т. 6. № 1. С. 3-9.
2. Барахнин В. Б., Ткачев Д. А. Кластеризация текстовых документов на основе составных ключевых термов // Вестник Новосибирского государственного университета. Серия: Информационные технологии. 2010. Т. 8. № 2. С. 5-14.
3. Баранов М. А. Модификация жадного алгоритма кластеризации // Прикладная информатика. 2013. № 2. С. 78-88.
4. Бериков В. С., Лобов Г. С. Современные тенденции в кластерном анализе // Всероссийский конкурсный отбор обзорно-аналитических статей по приоритетному направлению «Информационно-телекоммуникационные системы», 2008. С. 1-26.
5. Кормен Т., Лейзерсон Ч., Ривест Р., Штайн К. Алгоритмы: построение и анализ = Introduction to Algorithms / под ред. И. В. Красикова. 2-е изд. М.: Вильямс, 2005. — 1296 с.
6. Островский А. А. Реализация параллельного выполнения алгоритма FCM-кластеризации // Прикладная информатика. 2009. № 2. С. 101-106.
7. Nadizadeh Ali, Sahraeian Hashed, Zadeh Ali Sabze-vari, Homayouni Seyed Mahdi. Using greedy clustering method to solve capacitated location-routing problem // African Journal of Business Management Vol. 5 (21). 2011. Р. 79-85.
8. Chang D., Kantardzic M., Ouyang M. Hierarchical clustering with CUDA/GPU. In Proceedings of ISCA PDCCS. 2009. Р. 130-135.
9. Bohm C. et al. Density-based clustering using graphics processors // Proceedings of the 18th ACM conference on Information and knowledge management // ACM. 2009. С. 661-670.
10. Yue SH, Li P., Guo JD, Zhou SG. Using Greedy algorithm: DBSCAN revisited II // Journal of Zhejiang University SCIENCE. Vol. 5 (11). 2004. Р. 1405-1412.
M. Baranov, Graduate Student of National Research University «Higher School of Economics», Moscow, thenorthcat@gmail.com
parallel greedy clustering algorithm
The article is devoted to the development the parallel version of clustering algorithm based on the so-called greedy approach. Technology CUDA developed by NVIDIA is chosen for algorithm parallelization. The code and computing results for the similarity matrices of different sizes are presented.
Keywords: clustering, greedy algorithm, parallel computing, CUDA.
