Создание и практическое использование алгоритмов сильно зависит от возможности эффективно их реализовать. В лаборатории языковых инструментов JetBrains разрабатывают алгоритмы поиска путей в помеченных графах с дополнительными ограничениями. Эти алгоритмы достаточно естественно выражаются в терминах операций над булевыми матрицами, но в современных высокопроизводительных библиотеках линейной алгебры пока нет полного набора необходимых операций над булевым полукольцом. Поэтому мы решили их реализовать.
Всем привет! Меня зовут Мария Карпенко и в этом году я окончила корпоративную магистерскую программу JetBrains в Университете ИТМО по направлению Software Engineering. Я начала интересоваться программированием, обучаясь в СПбГЭУ на направлении «Математические методы в экономике», и на последнем курсе поступила в Computer Science Center (CSC). Именно в CSC я узнала о существовании корпоративной магистратуры и однозначно определилась с выбором программы. Сейчас я стажируюсь в «Яндексе» в команде Трекера.
В первый год магистратуры студенты активно приобретают новые знания, а на втором курсе начинается работа над дипломом. Я старалась выбрать для себя тему, которая позволит погрузиться в новые технологии и алгоритмы, с акцентом на практическую часть. В итоге я остановилась на теме «Реализация операций с разреженными булевыми матрицами на OpenCL». Работу выполняла под руководством Семёна Григорьева, исследователя лаборатории языковых инструментов JetBrains.
Операции над матрицами и графы
Многие алгоритмы для анализа графов можно выразить через операции с их матрицами смежности. Например, транзитивное замыкание, поиск в ширину, прямое произведение графов. Такие базовые матричные операции входят в стандарт GraphBLAS, разработанный для описания алгоритмов обработки графов в терминах линейной алгебры.
В большинстве практических задач графы оказываются сильно разреженными, что приводит к необходимости хранить матрицы и осуществлять операции с ними в сжатых форматах. В лаборатории языковых инструментов JetBrains разрабатывают алгоритмы анализа графов, сводимые к операциям над разреженными булевыми матрицами. Для подтверждения их практической применимости возникла необходимость в эффективной реализации следующих операций с разреженными булевыми матрицами:
матричное умножение;
матричное сложение;
транспонирование матрицы;
извлечение подматрицы;
редуцирование строк матрицы в вектор;
произведение Кронекера.
Специальные реализации этих операций для матриц, определенных над булевым полукольцом, позволяют, во-первых, экономить память: все форматы для хранения разреженных матриц содержат массив со значениями, который в случае булевой матрицы не нужен. Во-вторых, операции с элементами разреженных булевых матриц сводятся к определению отсутствия или наличия значения на конкретной позиции, что значительно быстрее операций с типами float или double.
Современные высокопроизводительные библиотеки линейной алгебры содержат лишь некоторые из необходимых операций, чаще всего умножение и сложение без специализации для булевых матриц. На момент начала работы библиотека cuBool, разрабатываемая в той же лаборатории языковых инструментов, содержала полный набор необходимых операций с реализацией на CUDA. На CPU отмечу библиотеку SuiteSparse, в которой полностью реализован стандарт GraphBLAS.
Моей задачей было создать полное решение на OpenCL, которое позволит задействовать видеокарты от AMD и другие устройства, поддерживающие стандарт. От меня не требовалось изобретать новые алгоритмы для операций с матрицами, необходимо было выбрать лучшие и обосновать свой выбор.
Знакомство с технологией и сжатыми форматами
Работа началась с погружения в высокопроизводительные вычисления, в чем мне помог курс от CSC «Вычисления на видеокартах». Помимо доступных и понятных объяснений принципов программирования на видеокартах, курс содержит ряд практических советов и приемов, связанных как с организацией кода, так и с поиском ошибок параллельного программирования. В рамках курса мы реализовывали многие необходимые промежуточные операции, которые пригодились мне в дальнейшем.
Также в начале работы мы определились с форматом хранения матриц. Так как нам заранее неизвестна структура матрицы, были выбраны наиболее общие и популярные форматы: COO и CSR (CSC). Именно с этими форматами работают большинство библиотек линейной алгебры. Так как в наших задачах матрицы могут быть сверхразреженными, мы сосредоточились на форматах, в которых матрица занимает O(nnz) памяти, где nnz — количество ненулевых элементов. Таким свойством обладают форматы COO и DCSR, полученный из CSR удалением информации о пустых строках.
Я реализовала операции в разных форматах, выбирая тот или иной формат исходя из особенностей работы с ним. Также были реализованы конвертации между форматами, которые осуществляются достаточно быстро. Некоторые операции по-разному можно реализовать как в COO, так и в DCSR, поэтому основную реализацию я выбирала по результатам замеров.
Выбранные форматы для операций:
COO:
транспонирование;
сложение;
произведение Кронекера.
DCSR:
умножение;
извлечение подматрицы;
редуцирование строк матрицы;
произведение Кронекера.
О реализациях операций
Первыми я реализовала сложение, транспонирование и произведение Кронекера в формате COO. Транспонирование в этом формате получается очень естественно — достаточно переставить массивы строковых и столбцовых индексов и отсортировать результат. Сложение в формате COO — это слияние двух отсортированных массивов с последующим удалением дубликатов.
Произведение Кронекера в COO также достаточно простое, но, как выяснилось позже, реализация через DCSR работает на порядок быстрее.
Для остальных операций однозначно был выбран формат DCSR, так как в них часто требуется индексация строк, которая в COO осуществляется медленнее.
Самая сложная из реализуемых операций — умножение. Алгоритмы умножения матриц развиваются уже несколько десятилетий, и я ознакомилась со многими работами в этой области. Среди всех выделяются работы W. Liu и Y. Nagasaka. Они подтвердили свою практическую применимость в библиотеках, а алгоритм Y. Nagasaka хорошо показал себя в cuBool.
Оба алгоритма опираются на формулировку умножения матриц, основанную на формировании строк матрицы-результата как линейной комбинации строк второй матрицы, то есть если A * B = C, то:
![C[i,:]=\sum_jA[i,j]*B[j,:]](https://habrastorage.org/getpro/habr/upload_files/69f/ee1/90d/69fee190db79bf2dcff3807e4bbe1351.svg)
Для булевых матриц задача сводится к формированию строк матрицы C посредством слияния некоторых строк матрицы B и последующего удаления дубликатов. Точно так же складываются матрицы в COO формате.
В алгоритме W. Liu слияние и удаление дубликатов осуществляются различными способами в зависимости от оценки размера строки. Среди них есть пирамидальная сортировка, битоническая сортировка с удалением дубликатов с помощью префиксной суммы, Merge Path. В алгоритме Y. Nagasaka используется один метод удаления дубликатов — хэш-таблицы с открытой адресацией, в которую обращаются несколько потоков одновременно.
Я реализовала оба варианта, и в моей реализации алгоритм Y. Nagasaka оказался лучше и по времени, и по памяти.
Отдельно хочется отметить, что алгоритм Y. Nagasaka значительно проще в поддержке и реализации: он состоит из однообразных ядер (программ для видеокарты), которые отличаются только распределением работы потоков, размером и размещением хэш-таблицы.
Извлечение подматрицы и редуцирование строк матрицы в вектор были реализованы последними в формате DCSR.
Эксперименты
В экспериментальной части мы сосредоточились на операциях, которые есть в большинстве библиотек линейной алгебры: сложение и умножение. Сложение мы тестировали на операции M + M2, а умножение — на M2. Для замеров были выбраны десять сильно разреженных матриц из коллекции SuiteSparse, но далее я приведу результаты по трем матрицам, которые можно считать наиболее показательными: две матрицы с равномерно распределенными ненулевыми элементами разных размеров (строки 1 и 3) и матрица с некоторыми плотными строками (строка 2).
Операции выполнялись на трех различных устройствах: две видеокарты (AMD, NVIDIA) и одна FPGA.
Характеристики устройств:
Вендор | NVIDIA | AMD | Intel FPGA |
Имя | GeForce GTX 1070 | Radeon Vega Frontier Edition | Arria 10 |
Глобальная память | 8 Гб | 16 Гб | 8 Гб |
Локальная память | 48 Кб | 64 Кб | 16 Кб |
Макс. размер блока | 1024 Б | 256 Б | 2 Гб |
Частота | 1746 MHz | 1600 MHz | 1000 MHz |
Число АЛУ | 1920 | 4096 | 427200 |
Скомпилировать решение на FPGA позволяет технология Intel FPGA SDK for OpenCL. За предоставленный сервер с устройством Arria10 благодарим компанию Selectel.
На FPGA мои реализации сгенерировались в неэффективную архитектуру, и этому есть несколько объяснений. Сами алгоритмы сильно заточены под устройство видеокарт: нагрузка распределяется с учетом организации работы потоков по ворпам и мультипроцессорам. Так как видеокарте не нужно перестраивать свою архитектуру под каждый бинарный файл, мы можем быстро скомпилировать множество ядер, суммарный объем оперативной памяти в которых ничем не ограничен, но в каждом отдельном ядре нужно учитывать ограничения мультипроцессора. С FPGA все иначе: ограничение в 16 КБ затрагивает всю программу, и мне не удалось, например, скомпилировать полный набор ядер для умножения. Пришлось уменьшать число ядер, что привело к неэффективному решению. Время выполнения обеих операций оказалась в сотни раз медленнее CPU-библиотеки SuiteSparse.
На операции сложения в моей реализации разница между результатами AMD и NVIDIA объясняется характеристиками видеокарт — AMD почти в 2 раза мощнее.
Сравнение выполнения операции M + M2 на различных устройствах:
№ | Число строк, млн | nnz(M), млн | nnz(M + M2), млн | NVIDIA, мс | AMD, мс | Intel FPGA, мс |
1 | 0,06 | 0,24 | 0,92 | 1,73 ± 3,97 | 1,71 ± 0,30 | 86,51 |
2 | 0,92 | 5,11 | 30,81 | 38,10 ± 0,53 | 19,34 ± 4,68 | 3077,81 |
3 | 2,22 | 4,88 | 13,63 | 15,51 ± 0,48 | 9,27 ± 2,81 | 1248,45 |
В умножении преимущества AMD становится заметно только на матрицах с более плотными строками:
Сравнение выполнения операции M2 на различных устройствах:
№ | Число строк, млн | nnz(M), млн | nnz(M2), млн | NVIDIA (1), мс | NVIDIA (2), мс | AMD (2), мс | Intel FPGA (2), мс |
1 | 0,06 | 0,24 | 0,71 | 2,71 ± 0,11 | 1,88 ± 0,34 | 2,78 ± 1,46 | 5590,79 |
2 | 0,92 | 5,11 | 29,71 | 184,06 ± 0,53 | 126,36 ± 0,34 | 85,97 ± 8,78 | 77601,9 |
3 | 2,22 | 4,88 | 8,76 | 57,27 ± 1,71 | 24,57 ± 0,11 | 29,66 ± 0,85 | - |
(1) — алгоритм W. Liu
(2) — алгоритм Y. Nagasaka
В таблицах ниже представлено сравнение с некоторыми библиотеками линейной алгебры на видеокарте NVIDIA. Как видно из замеров, сложение через перевод матриц в координатный формат реализуют все библиотеки, кроме cuBool. Отмечу, что по потреблению памяти cuBool серьезно выигрывает, так как в используемом алгоритме не выделяется дополнительная память на некоторые промежуточные вычисления.
Сравнение выполнения операции M + M2 среди библиотек линейной алгебры:
№ | nnz(M + M2), млн | cuBool, мс | clBool, мс | Cusp, мс | cuSPARSE, мс | SuiteSpars, мс |
1 | 0,92 | 1,12 ± 0,02 | 1,73 ± 3,97 | 1,54 ± 0,20 | 2,40 ± 0,04 | 3,91 ± 2,06 |
2 | 30,81 | 24,20 ± 0,90 | 38,10 ± 0,53 | 32,08 ± 1,34 | 88,48 ± 0,28 | 77,68 ± 1,66 |
3 | 13,63 | 30,65 ± 1,89 | 15,51 ± 0,48 | 14,75 ± 0,31 | 18,15 ± 0,03 | 50,29 ± 0,93 |
Несмотря на один алгоритм и похожую реализацию, на плотных матрицах cuBool работает заметно быстрее моей реализации clBool. Также cuBool и clBool требуют в среднем на треть меньше памяти, чем другие библиотеки.
Сравнение выполнения операции M2 среди библиотек линейной алгебры:
№ | nnz(M2), млн | cuBool, мс | clBool-hash, мс | clSPARSE, мс | cuSPARSE, мс | Cusp, мс | SuiteSparse, мс |
1 | 0,71 | 1,78 ± 0,06 | 1,88 ± 0,34 | 2,02 ± 0,34 | 2,04 ± 0,04 | 5,39 ± 0,12 | 7,92 ± 0,27 |
2 | 29,71 | 42,19 ± 1,26 | 126,36 ± 0,34 | 165,03 ± 8,30 | 4726,33 ± 0,52 | 246,232 ± 9,90 | 712,54 ± 20,07 |
3 | 8,76 | 35,30 ± 0,22 | 24,57 ± 0,11 | 38,51 ± 2,12 | 50,68 ± 0,14 | 51,28 ± 0,35 | 93,47 ± 2,09 |
Заключение
Результат моей работы используется в качестве OpenCL-бэкенда проекта spbla, который был создан для реализации, тестирования и изучения алгоритмов анализа графов, выразимых в терминах операций над булевыми матрицами.
Набор решаемых задач все еще серьезно ограничен по памяти: в ходе итераций алгоритма матрицы становятся плотнее и перестают помещаться в видеопамять. Естественное направление развития проекта — это распределенные алгоритмы для матричных операций и использование виртуальной памяти. Так как библиотеки с подобными решениями появятся нескоро, имеет смысл создавать решения в рамках лаборатории.
Хочется выразить благодарность лаборатории языковых инструментов за предоставление серверов с видеокартами для экспериментов. Отдельно выражаю благодарность научному руководителю Семёну Григорьеву за непрерывное и чуткое руководство, помощь в работе с текстом.
До 9 августа продолжается прием документов на корпоративную магистерскую программу JetBrains «Разработка программного обеспечения». В этом году набор ведется на 30 бюджетных и 5 платных мест. Подробнее о программе можно узнать на сайте или в Telegram-чате для абитуриентов.
