2025-11-20T19:52:15.672703

A GPU-resident Memory-Aware Algorithm for Accelerating Bidiagonalization of Banded Matrices

Ringoot, Alomairy, Edelman
The reduction of a banded matrix to a bidiagonal form is a crucial step in the Singular Value Decomposition (SVD), a cornerstone of scientific computing and AI. Despite being a highly parallel algorithm, it was previously believed to be unsuitable for GPU computation because it is memory bandwidth-bound. Recent developments in GPU hardware, including larger L1 memory per Streaming Multiprocessor/Compute Unit, have changed that. We present the first GPU algorithm for reducing a banded matrix to bidiagonal form as part of the NextLA.jl open-source software package. Our algorithm is based on previous CPU-based multicore parallel cache-efficient bulge chasing algorithms and adapted to optimize for GPU throughput. We leverage Julia Language's Array abstractions and KernelAbstractions to implement a single hardware- and data precision-agnostic function on NVIDIA, AMD, Intel, and Apple Metal GPUs for half, single, and double precision, and examine performance optimization across hardware architectures and data precision. We also develop a hardware-aware performance model and identify key hyperparameters, such as inner tilewidth and block concurrency, that govern optimal GPU execution for bandwidth-bound workloads. We demonstrate highly parallel bandwidth-bound algorithm on the GPU can outperform CPU-based implementations: the GPU algorithm outperforms multithreaded CPU High-Performance libraries PLASMA and SLATE as of matrix size 1024 x 1024 and by a factor over 100 for matrices of 32k x 32k. In addition, the performance of the algorithm increases linearly with matrix bandwidth size, making faster reduction of larger matrix bandwidths now also possible. With this work, we break memory bandwidth barriers, as well as matrix bandwidth barriers, resulting in orders-of-magnitude faster algorithms for the reduction of banded matrices to bidiagonal form on the GPU.
academic

GPU-резидентный алгоритм с учётом памяти для ускорения бидиагонализации ленточных матриц

Основная информация

  • ID статьи: 2510.12705
  • Название: A GPU-resident Memory-Aware Algorithm for Accelerating Bidiagonalization of Banded Matrices
  • Авторы: Evelyne Ringoot, Rabab Alomairy, Alan Edelman (MIT Computer Science & AI Laboratory)
  • Классификация: cs.DC (распределённые, параллельные и кластерные вычисления), cs.MS (математическое программное обеспечение)
  • Дата публикации: 14 октября 2025 г. (препринт arXiv)
  • Ссылка на статью: https://arxiv.org/abs/2510.12705

Аннотация

В данной работе предлагается первый GPU-резидентный алгоритм с учётом памяти для ускорения приведения ленточных матриц к бидиагональной форме, что является критическим этапом в разложении по сингулярным значениям (SVD). Несмотря на высокую параллелизуемость алгоритма, ранее считалось, что он непригоден для GPU-вычислений из-за ограничений пропускной способности памяти. Развитие GPU-оборудования, в частности увеличение объёма L1-памяти на каждый потоковый мультипроцессор/вычислительный блок, изменило эту ситуацию. Авторы адаптировали предыдущий алгоритм bulge-chasing, эффективный по кэшу для многоядерных CPU, и оптимизировали его для пропускной способности GPU. Алгоритм реализован в виде единственной функции, независимой от оборудования и точности данных, на GPU NVIDIA, AMD, Intel и Apple Metal, поддерживающей вычисления с половинной, одинарной и двойной точностью. Экспериментальные результаты показывают, что GPU-алгоритм превосходит многопоточные CPU-библиотеки высокого производства PLASMA и SLATE начиная с матриц размером 1024×1024, а на матрицах размером 32k×32k достигает ускорения более чем в 100 раз.

Научный контекст и мотивация

Определение проблемы

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

  1. Приведение плотной матрицы к ленточной форме
  2. Приведение ленточной матрицы к бидиагональной форме (bulge-chasing)
  3. Приведение бидиагональной матрицы к диагональной форме

Научная мотивация

Хотя GPU-реализации первого и третьего этапов широко изучены, второй этап остаётся недостаточно исследованным на современных GPU. Dongarra и соавторы в 2014 году отметили, что «ускорители плохо справляются с обработкой мелкозернистых вычислительных задач, ограниченных памятью (таких как bulge chasing), что ограничивает потенциальные выгоды от GPU-реализации второго этапа».

Технологические возможности

Недавние достижения в архитектуре GPU, в частности:

  • расширение объёма L1-кэша на каждый потоковый мультипроцессор
  • повышение пропускной способности внутрикристаллической памяти
  • более гибкие схемы доступа к памяти
  • улучшенное переиспользование кэша и расширенная иерархия памяти

Эти улучшения значительно изменили баланс между памятью и вычислениями, создав новые возможности для переработки алгоритмов, ограниченных памятью.

Основные вклады

  1. Первый GPU-резидентный алгоритм bulge-chasing с учётом памяти: предложен первый нативный GPU-алгоритм для приведения ленточной матрицы к бидиагональной форме, полностью использующий превосходные характеристики памяти новейшего поколения высокопроизводительных GPU, обеспечивающий ускорение в 10-100 раз по сравнению с оптимизированными CPU-реализациями.
  2. Стратегия кэш-эффективного bulge-chasing для матриц с большой полосой пропускания: благодаря новой кэш-эффективной стратегии постепенного сокращения полосы пропускания GPU-алгоритм может обрабатывать большие полосы пропускания, чем ранее, изменяя оптимальный компромисс полосы пропускания между первым и вторым этапами в крупномасштабном SVD.
  3. Открытая реализация, независимая от оборудования и точности данных: библиотека с открытым исходным кодом, реализованная на основе абстракций языка Julia, с единственной функцией, поддерживающей вычисления с одинарной, половинной и двойной точностью на GPU NVIDIA, AMD, Intel и Apple.

Описание методологии

Определение задачи

Дана ленточная матрица A ∈ ℝⁿˣⁿ с полосой пропускания BW. Целью является приведение её к бидиагональной матрице B посредством ортогональных преобразований таким образом, чтобы A = UBVᵀ, где U и V — ортогональные матрицы.

Архитектура алгоритма

Основной процесс алгоритма

Алгоритм 1: Приведение ленточной матрицы к бидиагональной форме 
с использованием векторов Хаусхолдера
Входные данные: полоса пропускания BW, внутренняя ширина плитки TW, размер матрицы n
1: for сокращение полосы пропускания i = (BW-1)/TW → 1 do
2:   целевая полоса пропускания TBW = 1 + i·TW  
3:   for параллельно: каждая строка R = 1→n do
4:     bulge строки k = R
5:     for каждый bulge строки: j = 0, j+=1 do
6:       if 3(R-1) < j и k ≤ n then
7:         вычислить вектор HH для строки k для исключения TW элементов 
         и применить к строкам ниже
8:         вычислить вектор HH для самого левого порождённого bulge столбца 
         и применить к столбцам справа
9:         обновить значение k
10:      end if
11:    end for
12:  end for
13: end for

Реализация GPU-ядра с учётом памяти

Алгоритм 2: GPU-ядро с учётом памяти, TPB потоков на блок
1: память потока: Ai (TW+1)
2: память блока: X (TW+1)  
3: все потоки в блоке совместно: X ← A[k,..]
4: синхронизация потоков
5: все потоки в блоке совместно: HH(X)
6: все потоки в блоке совместно: A[k,..]← X
7: синхронизация потоков
8: for l: 0→(CBW + TW)/TPB - 1 do
9:   поток i: вычислить строку r = k + l·CPB + i
10:  Ai ← A[r, ...]
11:  HH(X, Ai)  
12:  A[r, ...]← Ai
13: end for

Технические инновации

1. Стратегия разбиения по полосе пропускания

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

  • лучшее переиспользование кэша
  • улучшенную локальность памяти
  • снижение затрат на синхронизацию

2. Контролируемое планирование блоков

Введён параметр Max blocks, ограничивающий количество одновременных блоков потоков на каждый исполнительный блок GPU. Когда количество требуемых блоков bulge-chasing превышает лимит, применяется развёртывание цикла на уровне программного обеспечения, при котором один блок получает несколько задач и выполняет их последовательно в рамках одного запуска ядра.

3. Стратегия трёхцикловой развязки

Для обеспечения корректности и включения параллелизма алгоритм принудительно применяет трёхциклову развязку между последовательными сканированиями строк. То есть после завершения трёх bulge строк следующая строка может начать сканирование без перекрытия доступа к данным.

Экспериментальная установка

Аппаратная среда

  • CPU: Intel Xeon Platinum 8462Y+ 32 ядра, 64 потока
  • GPU: включая NVIDIA A100/H100/RTX4060, AMD MI250X/MI300X, Intel PVC 1100, Apple M1

Эталонные библиотеки

  • PLASMA (v25.5.27): пакет параллельной линейной алгебры
  • SLATE (v2024.10.29): масштабируемая библиотека линейной алгебры
  • CUBLAS geam: ядро добавления матриц в качестве эталона пропускной способности памяти

Метрики оценки

  • время выполнения и коэффициент ускорения
  • пропускная способность памяти (DRAM, L1, L2)
  • численная точность (относительная ошибка)
  • использование аппаратных ресурсов

Конфигурация тестирования

  • размер матрицы: 1024×1024 до 65k×65k
  • диапазон полосы пропускания: 32 до 512
  • точность данных: FP16, FP32, FP64

Экспериментальные результаты

Основные результаты производительности

Сравнение производительности GPU и CPU-библиотек

  • Малая полоса пропускания (32): GPU-реализация быстрее PLASMA в 100-1000 раз, быстрее SLATE в 10-300 раз
  • Большая полоса пропускания (512): GPU-версия быстрее SLATE в 2-9 раз, быстрее PLASMA в 0,9-6 раз
  • Эффект размера матрицы: производительность GPU растёт быстрее, чем CPU-реализация

Производительность на различном оборудовании

  • NVIDIA H100: базовая производительность
  • AMD MI300X: снижение производительности примерно в 1,5-2 раза
  • Intel PVC: снижение производительности примерно в 20 раз (несмотря на больший кэш)
  • Apple M1: хорошие результаты среди потребительских GPU

Результаты настройки гиперпараметров

Ключевые находки:

  1. Внутренняя ширина плитки является наиболее важным фактором производительности
    • оптимальное значение для FP32: 32 (соответствует 128-байтовой строке кэша)
    • оптимальное значение для FP64: 16 (соответствует 128-байтовой строке кэша)
  2. Важность параметров варьируется в зависимости от полосы пропускания:
    • полоса пропускания 32: максимальное количество блоков более критично
    • полоса пропускания 128: количество потоков на блок более критично

Анализ численной точности

  • FP64: ошибка близка к машинной точности
  • FP32: предсказуемый рост ошибки, зависящий от размера
  • FP16: сохраняет приемлемую точность на матрицах до 16k×16k

Влияние эволюции оборудования

  • MI300X показывает значительное улучшение по сравнению с MI250X (L1-кэш удвоен, добавлен унифицированный слой L2.5)
  • H100 продолжает превосходить A100 (L1-кэш увеличен на 33%, L2 на 25%)

Связанные работы

Развитие решателей SVD

  • CPU-реализации: PLASMA, SBR, FLAME, ELPA, LAPACK и другие, сосредоточенные на блочных преобразованиях Хаусхолдера и алгоритмах bulge-chasing
  • GPU-реализации: в основном сосредоточены на первом этапе (плотное к ленточному) или полном одноэтапном SVD, избегая нерегулярного доступа к памяти

Сравнение с решателями собственных значений

В отличие от решателей SVD, двухэтапные решатели для симметричных ленточных задач собственных значений на гибридных CPU-GPU системах имеют обширную литературу, однако исследования несимметричного приведения ленточной матрицы к бидиагональной форме отстают.

Алгоритм Bulge-chasing

Первоначально предложен как параллельный алгоритм, но исследования случая бидиагональной матрицы на CPU крайне редки. В сравнении, исследования bulge-chasing для симметричного приведения ленточной матрицы к трёхдиагональной форме и верхней форме Хессенберга более обширны.

Заключение и обсуждение

Основные выводы

  1. Преодоление барьера пропускной способности памяти: доказано, что алгоритмы линейной алгебры, ограниченные памятью, не только осуществимы на современных GPU, но и могут значительно превосходить производительность CPU
  2. Важность характеристик оборудования: отношение задержки L1/L2-кэша к размеру более критично, подчёркивая важность низколатентного быстрого доступа к памяти
  3. Принципы проектирования алгоритмов: оптимальная производительность достигается благодаря принципиальному проектированию алгоритмов, параметризации с учётом кэша и переносимой инфраструктуре компилятора

Ограничения

  1. Модель занятости GPU: алгоритм имеет значительно меньше блоков GPU на начальных этапах, чем доступных исполнительных блоков, требуя достаточно больших матриц для полного использования ресурсов GPU
  2. Переполнение регистров: высокая внутренняя ширина плитки может привести к переполнению регистров, влияя на производительность
  3. Зависимость от оборудования: значительные различия в производительности между различными архитектурами GPU требуют аппаратно-специфической настройки

Направления будущих исследований

  1. Оптимизация архитектуры: использование новых возможностей, таких как переполнение регистров общей памяти в CUDA 13.0
  2. Расширение алгоритма: исследование стратегий обработки больших полос пропускания
  3. Полный конвейер SVD: построение полностью GPU-резидентного конвейера SVD

Глубокая оценка

Преимущества

  1. Новаторский вклад: первая реализация GPU-резидентного алгоритма приведения ленточной матрицы к бидиагональной форме, заполняющая важный технологический пробел
  2. Высокая практическая ценность: значительное ускорение (10-1000 раз) обеспечивает важное практическое применение
  3. Технические инновации: умная стратегия разбиения по полосе пропускания и проектирование с учётом памяти соответствуют особенностям архитектуры GPU
  4. Высокая портативность: реализация на основе единого исходного кода на Julia поддерживает GPU различных производителей и множественную точность

Недостатки

  1. Недостаточный теоретический анализ: отсутствует теоретический анализ сложности алгоритма и доказательства сходимости
  2. Ограниченный диапазон тестирования: тестирование в основном на синтетических матрицах, отсутствует проверка на реальных сценариях приложений
  3. Сложность настройки параметров: требуется сложная настройка гиперпараметров для различного оборудования

Влияние

  1. Академическое значение: переопределяет границы производительности алгоритмов, ограниченных памятью, на GPU
  2. Практическая ценность: обеспечивает критическую технологическую поддержку для крупномасштабных научных вычислений и приложений AI
  3. Вклад в экосистему: открытая реализация способствует развитию кроссплатформенной GPU-экосистемы линейной алгебры

Применимые сценарии

  • крупномасштабные вычисления SVD матриц
  • научные вычисления и инженерное моделирование
  • снижение размерности и извлечение признаков в машинном обучении
  • GPU-приложения, требующие высокопроизводительной линейной алгебры

Библиография

Статья цитирует 86 соответствующих источников, охватывающих GPU-вычисления, алгоритмы линейной алгебры, экосистему языка Julia и другие важные работы в различных областях, обеспечивая прочную теоретическую основу для исследования.