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

Un Algorithme Conscient de la Mémoire Résident sur GPU pour Accélérer la Bidiagonalisation de Matrices Bandes

Informations Fondamentales

  • ID de l'article: 2510.12705
  • Titre: A GPU-resident Memory-Aware Algorithm for Accelerating Bidiagonalization of Banded Matrices
  • Auteurs: Evelyne Ringoot, Rabab Alomairy, Alan Edelman (Laboratoire d'Informatique et d'IA du MIT)
  • Classification: cs.DC (Informatique Distribuée, Parallèle et en Grappe), cs.MS (Logiciels Mathématiques)
  • Date de Publication: 14 octobre 2025 (prépublication arXiv)
  • Lien de l'article: https://arxiv.org/abs/2510.12705

Résumé

Cet article propose le premier algorithme résident sur GPU conscient de la mémoire pour accélérer la réduction de matrices bandes en matrices bidiagonales, une étape cruciale de la décomposition en valeurs singulières (SVD). Bien que cet algorithme soit hautement parallélisable, il était auparavant considéré comme inadapté au calcul sur GPU en raison de sa nature contrainte par la bande passante mémoire. Cette situation a changé avec l'évolution du matériel GPU, notamment l'augmentation de la mémoire L1 par multiprocesseur de flux/unité de calcul. Les auteurs s'appuient sur un algorithme précédent de poursuite de bombement (bulge-chasing) efficace en cache pour processeurs multicœurs CPU et l'optimisent pour le débit GPU. L'algorithme a été implémenté en tant que fonction unique indépendante du matériel et de la précision des données sur les GPU NVIDIA, AMD, Intel et Apple Metal, supportant les calculs en demi-précision, simple précision et double précision. Les expériences montrent que l'algorithme GPU surpasse les bibliothèques CPU haute performance multithread PLASMA et SLATE à partir d'une taille de matrice de 1024×1024, avec une amélioration de performance dépassant 100 fois sur des matrices de 32k×32k.

Contexte de Recherche et Motivation

Définition du Problème

La décomposition en valeurs singulières (SVD) est un outil numérique fondamental en calcul scientifique, apprentissage automatique et analyse de données, largement appliquée à l'analyse en composantes principales, l'indexation sémantique latente, l'approximation de faible rang et la complétion de matrices. La SVD sur matériel moderne à grande échelle suit généralement un processus en trois étapes :

  1. Réduction de matrice dense en matrice bande
  2. Réduction de matrice bande en matrice bidiagonale (bulge-chasing)
  3. Réduction de matrice bidiagonale en matrice diagonale

Motivation de la Recherche

Bien que les implémentations GPU des première et troisième étapes aient été largement étudiées, la deuxième étape reste insuffisamment explorée sur les GPU modernes. Dongarra et al. ont souligné en 2014 que « les accélérateurs fonctionnent mal avec les tâches de calcul à grain fin contraintes par la mémoire (comme le bulge chasing), limitant les gains potentiels des implémentations GPU pour la deuxième étape ».

Opportunités Techniques

Les progrès récents de l'architecture GPU, en particulier :

  • L'expansion de la taille du cache L1 par multiprocesseur de flux
  • Une bande passante intra-puce plus élevée
  • Des modèles d'accès mémoire plus flexibles
  • Une meilleure réutilisation du cache et une hiérarchie mémoire plus grande

Ces améliorations ont significativement modifié l'équilibre mémoire-calcul, créant de nouvelles opportunités pour la reconception d'algorithmes contraints par la mémoire.

Contributions Principales

  1. Premier algorithme bulge-chasing résident sur GPU conscient de la mémoire: Propose le premier algorithme natif GPU pour la réduction de matrices bandes en matrices bidiagonales, exploitant pleinement les propriétés mémoire supérieures des GPU haute performance de dernière génération, avec une amélioration de performance de 10 à 100 fois par rapport aux implémentations CPU optimisées.
  2. Stratégie bulge-chasing efficace en cache pour matrices à large bande: Grâce à une nouvelle stratégie efficace en cache de réduction progressive de la bande passante, l'algorithme GPU peut traiter des bandes passantes plus larges qu'auparavant, modifiant l'équilibre optimal de la bande passante entre la première et la deuxième étape dans les SVD de grande taille.
  3. Implémentation open-source indépendante du matériel et de la précision des données: Bibliothèque open-source basée sur les abstractions du langage Julia, avec une définition de fonction unique supportant les calculs en simple précision, demi-précision et double précision sur les GPU NVIDIA, AMD, Intel et Apple.

Détails de la Méthode

Définition de la Tâche

Étant donnée une matrice bande A ∈ ℝⁿˣⁿ de largeur de bande BW, l'objectif est de réduire celle-ci en matrice bidiagonale B par transformations orthogonales telles que A = UBVᵀ, où U et V sont des matrices orthogonales.

Architecture de l'Algorithme

Flux d'Algorithme Principal

Algorithme 1: Réduction de matrice bande en matrice bidiagonale 
             utilisant des vecteurs de Householder
Entrée: Largeur de bande BW, largeur de tuile interne TW, taille de matrice n
1: pour Réduction de bande passante i = (BW-1)/TW → 1 faire
2:   Bande passante cible TBW = 1 + i·TW  
3:   pour Parallèle: chaque ligne R = 1→n faire
4:     Bombement de ligne k = R
5:     pour Chaque bombement de ligne: j = 0, j+=1 faire
6:       si 3(R-1) < j et k ≤ n alors
7:         Calculer le vecteur HH pour la ligne k afin d'éliminer TW éléments 
          et appliquer aux lignes en dessous
8:         Calculer le vecteur HH pour le bombement de colonne généré le plus à gauche 
          et appliquer aux colonnes à droite
9:         Mettre à jour la valeur de k
10:      fin si
11:    fin pour
12:  fin pour
13: fin pour

Implémentation du Noyau GPU Conscient de la Mémoire

Algorithme 2: Noyau GPU conscient de la mémoire, TPB threads par bloc
1: Mémoire de thread: Ai (TW+1)
2: Mémoire de bloc: X (TW+1)  
3: Tous les threads du bloc collaborent: X ← A[k,..]
4: Synchroniser les threads
5: Tous les threads du bloc collaborent: HH(X)
6: Tous les threads du bloc collaborent: A[k,..]← X
7: Synchroniser les threads
8: pour l: 0→(CBW + TW)/TPB - 1 faire
9:   Thread i: Calculer la ligne r = k + l·CPB + i
10:  Ai ← A[r, ...]
11:  HH(X, Ai)  
12:  A[r, ...]← Ai
13: fin pour

Points d'Innovation Technique

1. Stratégie de Partitionnement de Bande Passante

Partitionner la matrice en tuiles de bande passante, exécuter la réduction par blocs de bande passante consécutifs plutôt que de réduire toute la bande passante à la fois. Cette stratégie de partitionnement réalise :

  • Une meilleure réutilisation du cache
  • Une localité mémoire améliorée
  • Une réduction des frais généraux de synchronisation

2. Ordonnancement de Bloc Contrôlé

Introduire un paramètre Max blocks limitant le nombre de blocs de threads concurrents exécutés par unité d'exécution GPU. Lorsque le nombre de blocs bulge-chasing requis par l'algorithme dépasse la limite, utiliser le déroulement de boucle au niveau logiciel, où un seul bloc se voit attribuer plusieurs tâches et les exécute séquentiellement dans un seul lancement de noyau.

3. Stratégie de Séparation Trois-Cycles

Pour assurer la correction et activer le parallélisme, l'algorithme impose une séparation trois-cycles entre les balayages de lignes consécutives. C'est-à-dire que chaque trois bombements de ligne complétés, la ligne suivante peut commencer le balayage sans chevaucher les accès aux données.

Configuration Expérimentale

Environnement Matériel

  • CPU: Intel Xeon Platinum 8462Y+ 32 cœurs 64 threads
  • GPU: Incluant NVIDIA A100/H100/RTX4060, AMD MI250X/MI300X, Intel PVC 1100, Apple M1

Références de Comparaison

  • PLASMA (v25.5.27): Paquetage d'Algèbre Linéaire Parallèle
  • SLATE (v2024.10.29): Bibliothèque d'Algèbre Linéaire Scalable
  • CUBLAS geam: Noyau d'addition de matrice comme référence de bande passante mémoire

Métriques d'Évaluation

  • Temps d'exécution et rapport d'accélération
  • Débit mémoire (DRAM, L1, L2)
  • Précision numérique (erreur relative)
  • Utilisation des ressources matérielles

Configuration de Test

  • Taille de matrice: 1024×1024 à 65k×65k
  • Plage de bande passante: 32 à 512
  • Précision des données: FP16, FP32, FP64

Résultats Expérimentaux

Résultats de Performance Principaux

Comparaison de Performance GPU vs Bibliothèques CPU

  • Petite bande passante (32): L'implémentation GPU est 100-1000 fois plus rapide que PLASMA, 10-300 fois plus rapide que SLATE
  • Grande bande passante (512): La version GPU est 2-9 fois plus rapide que SLATE, 0,9-6 fois plus rapide que PLASMA
  • Effet de la taille de matrice: La performance GPU croît plus rapidement avec la taille de matrice que les implémentations CPU

Performance Intermatériel

  • NVIDIA H100: Performance de référence
  • AMD MI300X: Environ 1,5-2 fois de dégradation de performance
  • Intel PVC: Environ 20 fois de dégradation de performance (malgré un cache plus grand)
  • Apple M1: Bonne performance parmi les GPU grand public

Résultats d'Optimisation des Hyperparamètres

Découvertes clés :

  1. La largeur de tuile interne est le facteur de performance le plus important
    • Valeur optimale FP32: 32 (correspondant à une ligne de cache de 128 octets)
    • Valeur optimale FP64: 16 (correspondant à une ligne de cache de 128 octets)
  2. L'importance des paramètres varie avec la bande passante:
    • Bande passante 32: Le nombre maximum de blocs est plus critique
    • Bande passante 128: Le nombre de threads par bloc est plus critique

Analyse de Précision Numérique

  • FP64: L'erreur est proche de la précision machine
  • FP32: Croissance d'erreur prévisible dépendante de la taille
  • FP16: Maintient une précision acceptable sur les matrices 16k×16k

Impact de l'Évolution Matérielle

  • MI300X montre une amélioration significative par rapport à MI250X (cache L1 doublé, ajout d'une couche L2.5 unifiée)
  • H100 continue de surpasser A100 (cache L1 augmenté de 33%, L2 augmenté de 25%)

Travaux Connexes

Développement des Solveurs SVD

  • Implémentations CPU: PLASMA, SBR, FLAME, ELPA, LAPACK, etc., se concentrant sur les transformations de Householder par blocs et les algorithmes bulge-chasing
  • Implémentations GPU: Principalement concentrées sur la première étape (dense vers bande) ou la SVD complète en une étape, évitant les accès mémoire irréguliers

Comparaison avec les Solveurs de Valeurs Propres

Contrairement aux solveurs SVD, les solveurs à deux étapes pour les problèmes de valeurs propres de matrices bandes symétriques ont de nombreux travaux sur les systèmes hybrides CPU-GPU, mais la recherche sur la réduction non-symétrique de bande vers bidiagonale est en retard.

Algorithme Bulge-Chasing

Initialement proposé comme algorithme parallèle, mais la recherche sur le cas bidiagonal sur CPU est extrêmement rare. En comparaison, la recherche sur le bulge-chasing pour la réduction de bande symétrique vers tridiagonale et la forme supérieure de Hessenberg est plus extensive.

Conclusion et Discussion

Conclusions Principales

  1. Surmonter la Barrière de Bande Passante Mémoire: Démontre que les algorithmes d'algèbre linéaire contraints par la mémoire ne sont non seulement réalisables sur les GPU modernes, mais peuvent aussi surpasser significativement la performance CPU
  2. Importance des Caractéristiques Matérielles: La latence du cache L1/L2 est plus critique que la taille, soulignant l'importance de l'accès mémoire rapide à faible latence
  3. Principes de Conception d'Algorithmes: La performance optimale provient de la conception d'algorithmes principielle, de la paramétrisation consciente du cache et de l'infrastructure de compilateur portable

Limitations

  1. Modèle d'Occupation GPU: L'algorithme a un nombre de blocs GPU significativement inférieur aux unités d'exécution disponibles dans les étapes initiales, nécessitant des matrices suffisamment grandes pour utiliser pleinement les ressources GPU
  2. Débordement de Registres: Une largeur de tuile interne élevée peut entraîner un débordement de registres, affectant la performance
  3. Dépendance Matérielle: Les variations de performance entre architectures GPU différentes sont significatives, nécessitant une optimisation spécifique au matériel

Directions Futures

  1. Optimisation d'Architecture: Exploiter les nouvelles fonctionnalités comme le débordement de registres de mémoire partagée dans CUDA 13.0
  2. Extension d'Algorithme: Explorer les stratégies de traitement pour des bandes passantes plus larges
  3. Pipeline SVD Complet: Construire un pipeline SVD complètement résident sur GPU

Évaluation Approfondie

Points Forts

  1. Contribution Pionnière: Première implémentation d'un algorithme résident sur GPU pour la réduction bande vers bidiagonale, comblant une lacune technologique importante
  2. Valeur Pratique Élevée: L'amélioration de performance significative (10-1000 fois) la rend d'une grande valeur applicative
  3. Innovation Technique: La conception astucieuse du partitionnement de bande passante et de la conscience mémoire s'adapte bien aux caractéristiques de l'architecture GPU
  4. Portabilité Forte: L'implémentation à source unique basée sur Julia supporte les GPU multi-fournisseurs et multi-précision

Insuffisances

  1. Analyse Théorique Insuffisante: Manque d'analyse de la complexité algorithmique et de preuves de convergence
  2. Portée de Test Limitée: Principalement testée sur des matrices synthétiques, manquant de validation dans des scénarios d'application réelle
  3. Complexité d'Optimisation des Paramètres: Nécessite une optimisation complexe des hyperparamètres pour différents matériels

Impact

  1. Signification Académique: Redéfinit les limites de performance des algorithmes contraints par la mémoire sur GPU
  2. Valeur Pratique: Fournit un support technologique clé pour le calcul matriciel à grande échelle et les applications d'IA
  3. Contribution à l'Écosystème: L'implémentation open-source promeut le développement de l'écosystème d'algèbre linéaire GPU multi-plateforme

Scénarios Applicables

  • Calcul SVD de matrices à grande échelle
  • Calcul scientifique et simulation d'ingénierie
  • Réduction de dimensionnalité et extraction de caractéristiques en apprentissage automatique
  • Applications GPU nécessitant une algèbre linéaire haute performance

Références

L'article cite 86 références connexes, couvrant plusieurs domaines incluant le calcul GPU, les algorithmes d'algèbre linéaire, l'écosystème du langage Julia, etc., fournissant une base théorique solide pour la recherche.