Параллельные вычисления – матричное умножение
Сравнение производительности методов параллельных вычислений: матричное произведение
Введение
Матричное произведение является одной из центральных операций в численных методах, машинном обучении, обработке сигналов и физическом моделировании. В отличие от поэлементных операций над векторами, здесь вычислительная сложность растёт как , а узким местом становится не только пропускная способность памяти, но и эффективность использования кэша, векторных регистров и оптимизированных библиотек линейной алгебры.
Блокировка по кэшу – это приём организации вычислений, при котором большие матрицы разбиваются на небольшие блоки такого размера, чтобы они целиком помещались в кэш-память процессора. Вместо того чтобы многократно обращаться к медленной оперативной памяти, программа переиспользует блоки, уже загруженные в быстрый кэш, и таким образом сокращает время доступа к данным.
Переиспользование данных – это свойство алгоритма, при котором одни и те же элементы входных массивов используются многократно в ходе вычислений, пока они ещё находятся в кэше или регистрах. Чем выше переиспользование, тем меньше доля обращений к оперативной памяти и тем ближе производительность к пиковой вычислительной способности процессора.
В данной работе сравнивается производительность различных способов вычисления произведения двух квадратных матриц размера : от последовательного трёхциклового варианта до автоматической векторизации, тензорной нотации и оптимизированного вычисления.
Подключение библиотек
Присоединим необходимые библиотеки:
- BenchmarkTools — эталонный замер времени выполнения и аллокаций (
@btime); - LinearAlgebra — стандартный
mul!, вызывающий оптимизированный базовые подпрограммы линейной алгебры (BLAS); - Octavian — многопоточное BLAS-подобное умножение матриц;
- LoopVectorization — автоматическая векторизация циклов (
@turbo); - Tullio — тензорная нотация Эйнштейна с автоматическим распараллеливанием и векторизацией.
for pkg in (["BenchmarkTools", "LinearAlgebra", "Octavian", "LoopVectorization", "Tullio"])
try
eval(Meta.parse("using $pkg"))
catch
Pkg.add(pkg)
eval(Meta.parse("using $pkg"))
end
end
Постановка задачи
В качестве эталонной задачи выберем произведение двух квадратных матриц и размера .
Каждый элемент результата является скалярным произведением -й строки на -й столбец :
что в матричной форме записывается как .
Исходные данные
Создадим две случайные квадратные матрицы и размера и матрицу для результата. Размер выберем так, чтобы время последовательной реализации оставалось измеримым, но не слишком большим для интерактивной работы.
N = 1000
A = randn(N, N)
B = randn(N, N)
C = similar(A);
Сравнение методов
Полностью последовательное вычисление
Последовательная реализация через три вложенных цикла, обёрнутая в функцию для стабильной типизации JIT-компилятором. Порядок циклов j → i → k обеспечивает последовательный доступ к столбцам B, а аккумулятор s позволяет не создавать в C на каждой итерации внутреннего цикла.
function последовательное(C, A, B)
@inbounds for j in axes(C, 2)
for i in axes(C, 1)
s = zero(eltype(C))
for k in axes(A, 2)
s += A[i, k] * B[k, j]
end
C[i, j] = s
end
end
return C
end
@btime последовательное(C, A, B);
Время выполнения составило 3.698 с при нулевых аллокациях. Это базовая точка отсчёта, относительно которой будем считать ускорение остальных методов.
LoopVectorization
Тот же трёхцикловый алгоритм, но с макросом @turbo, который автоматически применяет SIMD-векторизацию (способ выполнения вычислений, при котором одна процессорная инструкция применяется сразу к нескольким элементам данных одновременно), разворачивание циклов и перестановку инструкций.
function turbo_matmul(C, A, B)
@turbo for j in axes(C, 2)
for i in axes(C, 1)
C[i, j] = zero(eltype(C))
for k in axes(A, 2)
C[i, j] += A[i, k] * B[k, j]
end
end
end
return C
end
@btime turbo_matmul(C, A, B);
Время составило 64.14 мс при нулевых аллокациях, что в ≈57.7 раз быстрее последовательного варианта. Однако накопление в C[i, j] внутри внутреннего цикла создаёт лишние обращения к памяти по сравнению с аккумулятором в регистре.
Octavian
Octavian это реализация матричного умножения, использующая многопоточность и блокировку по кэшу. Она не зависит от внешнего BLAS и в ряде случаев показывает результаты, конкурентные OpenBLAS, особенно на небольших и средних размерах матриц.
@btime Octavian.matmul!(C, A, B);
Время – 34.55 мс при нулевых аллокациях, ускорение в ≈107 раз относительно последовательного варианта.
Tullio
Нотация Эйнштейна: @tullio C[i, j] = A[i, k] * B[k, j]. По умолчанию Tullio применяет многопоточность и векторизацию, а также пытается использовать блокировку по кэшу и переставлять порядок циклов. Для матричного произведения это даёт результат, заметно превосходящий чистый @turbo.
@btime @tullio C[i, j] = A[i, k] * B[k, j];
Время — 12.99 мс при минимальных аллокациях (2.22 KiB на служебные структуры). Ускорение в ≈284.6 раз относительно последовательного варианта.
LinearAlgebra.mul!
Функция mul! из LinearAlgebra для Matrix{Float64} вызывает оптимизированную BLAS-реализацию. Это многопоточная, векторизованная и блокированная по кэшу операция.
@btime LinearAlgebra.mul!(C, A, B);
Время – 8.841 мс при нулевых аллокациях, ускорение в ≈418.3 раз относительно последовательного варианта. Это лучший результат среди всех измеренных методов на CPU.
Итоговая таблица
|
Место |
Метод |
Время |
Ускорение |
|
1 |
|
8.841 мс |
≈418× |
|
2 |
Tullio |
12.99 мс |
≈285× |
|
3 |
Octavian |
34.55 мс |
≈107× |
|
4 |
LoopVectorization |
64.14 мс |
≈58× |
|
5 |
Последовательное |
3698 мс |
1× |
Разброс производительности между последовательной реализацией и оптимизированным BLAS составляет около 418 раз.
Заключение
Проведённые измерения показывают, что для матричного произведения выбор метода даёт разброс производительности почти в три порядка: от 3.7 секунд последовательного трёхциклового варианта до 8.8 миллисекунд стандартной оптимизированной функции LinearAlgebra.mul!. При этом ключевой вывод состоит в том, что здесь простого распараллеливания или автовекторизации недостаточно. Последовательный @turbo ускоряет выполнение примерно в 58 раз, но всё равно проигрывает LinearAlgebra.mul! почти на порядок, потому что не выполняет блокировку по кэшу и переиспользование данных.
Промежуточное положение занимают Octavian и Tullio. Octavian показывает результат, близкий к LinearAlgebra.mul!, но на данной конфигурации немного ему уступает; Tullio, благодаря перестановке циклов и многопоточности, оказывается быстрее Octavian и лишь в 1.5 раза медленнее mul!. Это делает Tullio разумным выбором, если по каким-то причинам нельзя использовать стандартную библиотеку линейной алгебры, но нужна производительность, близкая к оптимальной.
Практический вывод: для матричного произведения в большинстве случаев наиболее эффективным инструментов будет являться LinearAlgebra.mul!, который автоматически вызывают высокооптимизированную библиотеку линейной алгебры.