Сравнение производительности методов параллельных вычислений: матричное произведение

Введение

Матричное произведение является одной из центральных операций в численных методах, машинном обучении, обработке сигналов и физическом моделировании. В отличие от поэлементных операций над векторами, здесь вычислительная сложность растёт как , а узким местом становится не только пропускная способность памяти, но и эффективность использования кэша, векторных регистров и оптимизированных библиотек линейной алгебры.

Блокировка по кэшу – это приём организации вычислений, при котором большие матрицы разбиваются на небольшие блоки такого размера, чтобы они целиком помещались в кэш-память процессора. Вместо того чтобы многократно обращаться к медленной оперативной памяти, программа переиспользует блоки, уже загруженные в быстрый кэш, и таким образом сокращает время доступа к данным.

Переиспользование данных – это свойство алгоритма, при котором одни и те же элементы входных массивов используются многократно в ходе вычислений, пока они ещё находятся в кэше или регистрах. Чем выше переиспользование, тем меньше доля обращений к оперативной памяти и тем ближе производительность к пиковой вычислительной способности процессора.

В данной работе сравнивается производительность различных способов вычисления произведения двух квадратных матриц размера : от последовательного трёхциклового варианта до автоматической векторизации, тензорной нотации и оптимизированного вычисления.

Подключение библиотек

Присоединим необходимые библиотеки:

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 s (0 allocations: 0 bytes)

Время выполнения составило 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.135 ms (0 allocations: 0 bytes)

Время составило 64.14 мс при нулевых аллокациях, что в ≈57.7 раз быстрее последовательного варианта. Однако накопление в C[i, j] внутри внутреннего цикла создаёт лишние обращения к памяти по сравнению с аккумулятором в регистре.

Octavian

Octavian это реализация матричного умножения, использующая многопоточность и блокировку по кэшу. Она не зависит от внешнего BLAS и в ряде случаев показывает результаты, конкурентные OpenBLAS, особенно на небольших и средних размерах матриц.

@btime Octavian.matmul!(C, A, B);

  34.552 ms (0 allocations: 0 bytes)

Время – 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.993 ms (59 allocations: 2.22 KiB)

Время — 12.99 мс при минимальных аллокациях (2.22 KiB на служебные структуры). Ускорение в ≈284.6 раз относительно последовательного варианта.

LinearAlgebra.mul!

Функция mul! из LinearAlgebra для Matrix{Float64} вызывает оптимизированную BLAS-реализацию. Это многопоточная, векторизованная и блокированная по кэшу операция.

@btime LinearAlgebra.mul!(C, A, B);

  8.841 ms (0 allocations: 0 bytes)

Время – 8.841 мс при нулевых аллокациях, ускорение в ≈418.3 раз относительно последовательного варианта. Это лучший результат среди всех измеренных методов на CPU.

Итоговая таблица

Место	Метод	Время	Ускорение
1	`LinearAlgebra.mul!`	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!, который автоматически вызывают высокооптимизированную библиотеку линейной алгебры.