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

Ein GPU-residenter speichergerichteter Algorithmus zur Beschleunigung der Bidiagonalisierung von Bandmatrizen

Grundinformationen

  • Papier-ID: 2510.12705
  • Titel: A GPU-resident Memory-Aware Algorithm for Accelerating Bidiagonalization of Banded Matrices
  • Autoren: Evelyne Ringoot, Rabab Alomairy, Alan Edelman (MIT Computer Science & AI Laboratory)
  • Klassifizierung: cs.DC (Verteiltes, paralleles und Cluster-Computing), cs.MS (Mathematische Software)
  • Veröffentlichungsdatum: 14. Oktober 2025 (arXiv-Preprint)
  • Papierlink: https://arxiv.org/abs/2510.12705

Zusammenfassung

Dieses Papier präsentiert den ersten GPU-residenten speichergerichteten Algorithmus zur Beschleunigung der Reduktion von Bandmatrizen zu Bidiagonalmatrizen, einem kritischen Schritt bei der Singulärwertzerlegung (SVD). Obwohl dieser Algorithmus hochgradig parallelisierbar ist, wurde er aufgrund seiner speicherbandbreitenbeschränkten Natur zuvor als ungeeignet für GPU-Berechnungen angesehen. Mit der Entwicklung der GPU-Hardware, insbesondere größerem L1-Speicher pro Stream-Multiprozessor/Recheneinheit, hat sich dies geändert. Die Autoren bauen auf dem vorherigen CPU-Multi-Core-Parallel-Cache-effizienten Bulge-Chasing-Algorithmus auf und optimieren ihn für GPU-Durchsatz. Der Algorithmus wurde als einzelne hardwareunabhängige und datenpräzisions-unabhängige Funktion auf NVIDIA-, AMD-, Intel- und Apple Metal-GPUs implementiert und unterstützt Halb-, Einfach- und Doppelgenauigkeitsberechnungen. Experimente zeigen, dass der GPU-Algorithmus bei Matrixgrößen von 1024×1024 die Multi-Thread-CPU-Hochleistungsbibliotheken PLASMA und SLATE übertrifft und bei 32k×32k-Matrizen eine Leistungssteigerung von über 100-fach erreicht.

Forschungshintergrund und Motivation

Problemdefinition

Die Singulärwertzerlegung (SVD) ist ein grundlegendes numerisches Werkzeug in wissenschaftlichen Berechnungen, maschinellem Lernen und Datenanalyse mit breiter Anwendung in Hauptkomponentenanalyse, latenter semantischer Indizierung, Niedrigrang-Approximation und Matrixvervollständigung. Moderne SVD auf großskaliger Hardware folgt typischerweise einem dreistufigen Prozess:

  1. Reduktion von dichter Matrix zu Bandmatrix
  2. Reduktion von Bandmatrix zu Bidiagonalmatrix (Bulge-Chasing)
  3. Reduktion von Bidiagonalmatrix zu Diagonalmatrix

Forschungsmotivation

Während GPU-Implementierungen der ersten und dritten Stufe bereits umfassend erforscht wurden, bleibt die zweite Stufe auf modernen GPUs noch unzureichend erforscht. Dongarra et al. wiesen 2014 darauf hin, dass „Beschleuniger bei der Verarbeitung speicherbeschränkter feinkorniger Berechnungsaufgaben (wie Bulge-Chasing) schlecht abschneiden und das Potenzial von GPU-Implementierungen der zweiten Stufe begrenzen".

Technische Chancen

Jüngste Fortschritte in der GPU-Architektur, insbesondere:

  • Erweiterung der L1-Cache-Größe pro Stream-Multiprozessor
  • Höhere On-Chip-Bandbreite
  • Flexiblere Speicherzugriffsmuster
  • Bessere Cache-Wiederverwendung und größere Speicherhierarchie

Diese Verbesserungen verändern das Speicher-Berechnungs-Gleichgewicht erheblich und schaffen neue Chancen für die Neugestaltung speicherbeschränkter Algorithmen.

Kernbeiträge

  1. Erster GPU-residenter speichergerichteter Bulge-Chasing-Algorithmus: Präsentiert den ersten GPU-nativen Algorithmus für die Reduktion von Bandmatrix zu Bidiagonalmatrix, der die überlegenen Speichereigenschaften der neuesten Generation hochleistungsfähiger GPUs vollständig nutzt und eine Leistungssteigerung von 10-100-fach gegenüber optimierten CPU-Implementierungen erreicht.
  2. Cache-effiziente Bulge-Chasing-Strategie für große Bandbreitematrizen: Durch eine neue Cache-effiziente Strategie der schrittweisen Bandbreitenreduktion kann der GPU-Algorithmus größere Bandbreiten als zuvor verarbeiten und verändert das optimale Bandbreitenverhältnis zwischen der ersten und zweiten Stufe in großen SVDs.
  3. Open-Source-Hardware-unabhängige und datenpräzisions-unabhängige Implementierung: Eine auf Julia-Sprachabstraktion basierende Open-Source-Bibliothek mit einer einzelnen Funktionsdefinition, die Einfach-, Halb- und Doppelgenauigkeitsberechnungen auf NVIDIA-, AMD-, Intel- und Apple-GPUs unterstützt.

Methodische Details

Aufgabendefinition

Gegeben eine Bandmatrix A ∈ ℝⁿˣⁿ mit Bandbreite BW besteht das Ziel darin, diese durch orthogonale Transformationen zu einer Bidiagonalmatrix B zu reduzieren, so dass A = UBVᵀ, wobei U und V orthogonale Matrizen sind.

Algorithmus-Architektur

Kern-Algorithmus-Ablauf

Algorithmus 1: Reduktion von Bandmatrix zu Bidiagonalmatrix mit Householder-Vektoren
Eingabe: Bandbreite BW, interne Kachelbreite TW, Matrixgröße n
1: für Bandbreitenreduktion i = (BW-1)/TW → 1 do
2:   Zielband TBW = 1 + i·TW  
3:   für parallel: jede Zeile R = 1→n do
4:     Zeilen-Bulge k = R
5:     für jeden Zeilen-Bulge: j = 0, j+=1 do
6:       wenn 3(R-1) < j und k ≤ n dann
7:         Berechne HH-Vektor für Zeile k, um TW Elemente zu eliminieren und auf untere Zeilen anzuwenden
8:         Berechne HH-Vektor für linksten erzeugten Spalten-Bulge und wende auf rechte Spalten an
9:         Aktualisiere k-Wert
10:      end wenn
11:    end für
12:  end für
13: end für

Speichergerichtete GPU-Kernel-Implementierung

Algorithmus 2: Speichergerichteter GPU-Kernel mit TPB Threads pro Block
1: Thread-Speicher: Ai (TW+1)
2: Block-Speicher: X (TW+1)  
3: Alle Threads im Block zusammenarbeiten: X ← A[k,..]
4: Thread-Synchronisation
5: Alle Threads im Block zusammenarbeiten: HH(X)
6: Alle Threads im Block zusammenarbeiten: A[k,..]← X
7: Thread-Synchronisation
8: für l: 0→(CBW + TW)/TPB - 1 do
9:   Thread i: Berechne Zeile r = k + l·CPB + i
10:  Ai ← A[r, ...]
11:  HH(X, Ai)  
12:  A[r, ...]← Ai
13: end für

Technische Innovationen

1. Bandbreitenkachelung-Strategie

Teilt die Matrix in Bandbreitenkacheln auf und führt die Reduktion nach aufeinanderfolgenden Bandbreitblöcken durch, anstatt die gesamte Bandbreite auf einmal zu reduzieren. Diese Kachelungsstrategie ermöglicht:

  • Bessere Cache-Wiederverwendung
  • Verbesserte Speicherlokalisierung
  • Reduzierte Synchronisationskosten

2. Gesteuerte Block-Planung

Führt einen Max-Blocks-Parameter ein, der die Anzahl gleichzeitiger Thread-Blöcke pro GPU-Ausführungseinheit begrenzt. Wenn die vom Algorithmus erforderliche Anzahl von Bulge-Chasing-Blöcken das Limit überschreitet, wird Software-Level-Loop-Unrolling verwendet, wobei einem einzelnen Block mehrere Aufgaben zugewiesen werden und diese sequenziell innerhalb eines einzelnen Kernel-Starts ausgeführt werden.

3. Drei-Zyklus-Trennungsstrategie

Um Korrektheit zu gewährleisten und Parallelität zu ermöglichen, erzwingt der Algorithmus eine Drei-Zyklus-Trennung zwischen aufeinanderfolgenden Zeilenscans. Das heißt, nach Abschluss von drei Zeilen-Bulges kann die nächste Zeile mit dem Scannen beginnen, ohne dass sich Datenzugriffe überlappen.

Experimentelle Einrichtung

Hardware-Umgebung

  • CPU: Intel Xeon Platinum 8462Y+ 32-Kern 64-Thread
  • GPU: Einschließlich NVIDIA A100/H100/RTX4060, AMD MI250X/MI300X, Intel PVC 1100, Apple M1

Vergleichsmaßstäbe

  • PLASMA (v25.5.27): Parallel Linear Algebra Software Package
  • SLATE (v2024.10.29): Scalable Linear Algebra Library
  • CUBLAS geam: Matrix-Addition-Kernel als Speicherbandbreitenreferenz

Bewertungsmetriken

  • Laufzeit und Beschleunigungsverhältnis
  • Speicherdurchsatz (DRAM, L1, L2)
  • Numerische Genauigkeit (relativer Fehler)
  • Hardware-Ressourcenauslastung

Test-Konfiguration

  • Matrixgröße: 1024×1024 bis 65k×65k
  • Bandbreitenbereich: 32 bis 512
  • Datenpräzision: FP16, FP32, FP64

Experimentelle Ergebnisse

Hauptleistungsergebnisse

GPU vs. CPU-Bibliotheken-Leistungsvergleich

  • Kleine Bandbreite (32): GPU-Implementierung 100-1000-fach schneller als PLASMA, 10-300-fach schneller als SLATE
  • Große Bandbreite (512): GPU-Version 2-9-fach schneller als SLATE, 0,9-6-fach schneller als PLASMA
  • Matrixgrößen-Effekt: GPU-Leistung wächst schneller mit Matrixgröße als CPU-Implementierungen

Leistung über Hardware hinweg

  • NVIDIA H100: Basis-Leistung
  • AMD MI300X: Etwa 1,5-2-fache Leistungsabnahme
  • Intel PVC: Etwa 20-fache Leistungsabnahme (trotz größerem Cache)
  • Apple M1: Gute Leistung unter Consumer-GPUs

Hyperparameter-Optimierungsergebnisse

Wichtigste Erkenntnisse:

  1. Interne Kachelbreite ist der wichtigste Leistungsfaktor
    • FP32 optimaler Wert: 32 (entspricht 128-Byte-Cache-Zeile)
    • FP64 optimaler Wert: 16 (entspricht 128-Byte-Cache-Zeile)
  2. Parameterimportanz variiert mit Bandbreite:
    • Bandbreite 32: Maximale Block-Anzahl kritischer
    • Bandbreite 128: Threads pro Block kritischer

Numerische Genauigkeitsanalyse

  • FP64: Fehler nahe Maschinengenauigkeit
  • FP32: Vorhersehbarer größenabhängiger Fehleranstieg
  • FP16: Behält akzeptable Genauigkeit bei 16k×16k-Matrizen

Auswirkungen der Hardware-Entwicklung

  • MI300X zeigt signifikante Verbesserung gegenüber MI250X (L1-Cache verdoppelt, einheitliche L2.5-Schicht hinzugefügt)
  • H100 übertrifft kontinuierlich A100 (L1-Cache +33%, L2 +25%)

Verwandte Arbeiten

SVD-Solver-Entwicklung

  • CPU-Implementierungen: PLASMA, SBR, FLAME, ELPA, LAPACK konzentrieren sich auf Block-Householder-Transformationen und Bulge-Chasing-Algorithmen
  • GPU-Implementierungen: Konzentrieren sich hauptsächlich auf die erste Stufe (dicht zu Band) oder vollständige einstufige SVD, vermeiden unregelmäßige Speicherzugriffe

Eigenwert-Solver-Vergleich

Im Gegensatz zu SVD-Solvern gibt es umfangreiche Arbeiten zu zweistufigen Solvern für symmetrische Band-Eigenwertprobleme auf hybriden CPU-GPU-Systemen, aber die Forschung zu nicht-symmetrischen Band-zu-Bidiagonal-Reduktionen hinkt hinterher.

Bulge-Chasing-Algorithmen

Ursprünglich als paralleler Algorithmus vorgeschlagen, aber die Forschung zum Bidiagonal-Fall auf CPUs ist äußerst selten. Im Vergleich dazu ist die Forschung zu symmetrischen Band-zu-Tridiagonal- und oberen Hessenberg-Form-Bulge-Chasing umfangreicher.

Schlussfolgerungen und Diskussion

Hauptschlussfolgerungen

  1. Durchbruch bei Speicherbandbreitenbeschränkungen: Beweist, dass speicherbeschränkte lineare Algebra-Algorithmen auf modernen GPUs nicht nur machbar, sondern auch erheblich CPU-Leistung übertreffen können
  2. Bedeutung von Hardware-Eigenschaften: L1/L2-Cache-Latenz ist wichtiger als Größe, unterstreicht die Bedeutung von niedriger Latenz und schnellem Speicherzugriff
  3. Algorithmus-Design-Prinzipien: Optimale Leistung ergibt sich aus prinzipiellem Algorithmus-Design, Cache-bewusster Parametrisierung und tragbarer Compiler-Infrastruktur

Einschränkungen

  1. GPU-Auslastungsmodell: Der Algorithmus hat in frühen Phasen deutlich weniger GPU-Blöcke als verfügbare Ausführungseinheiten, benötigt ausreichend große Matrizen zur vollständigen GPU-Auslastung
  2. Register-Spillover: Hohe interne Kachelbreite kann zu Register-Spillover führen und Leistung beeinträchtigen
  3. Hardware-Abhängigkeit: Signifikante Leistungsunterschiede zwischen GPU-Architekturen erfordern Hardware-spezifische Optimierung

Zukünftige Richtungen

  1. Architektur-Optimierung: Nutzung neuer Funktionen wie CUDA 13.0 mit Shared-Memory-Register-Spillover
  2. Algorithmus-Erweiterung: Erkundung von Verarbeitungsstrategien für größere Bandbreiten
  3. Vollständige SVD-Pipeline: Konstruktion einer vollständig GPU-residenten SVD-Pipeline

Tiefgreifende Bewertung

Stärken

  1. Bahnbrechender Beitrag: Erste GPU-residente Implementierung des Band-zu-Bidiagonal-Reduktions-Algorithmus, füllt wichtige technische Lücke
  2. Hoher praktischer Wert: Signifikante Leistungssteigerung (10-1000-fach) hat wichtigen Anwendungswert
  3. Technische Innovation: Geschickte Bandbreitenkachelung und speichergerichtetes Design passen zu GPU-Architektur-Charakteristiken
  4. Starke Portabilität: Julia-basierte Single-Source-Implementierung unterstützt Multi-Vendor-GPUs und Multi-Präzision

Mängel

  1. Unzureichende theoretische Analyse: Fehlt an Algorithmus-Komplexitätsanalyse und Konvergenzbeweis
  2. Begrenzte Test-Reichweite: Hauptsächlich auf synthetischen Matrizen getestet, fehlt Validierung in realen Anwendungsszenarien
  3. Komplexe Parameteroptimierung: Erfordert komplexe Hyperparameter-Optimierung für verschiedene Hardware

Einfluss

  1. Akademische Bedeutung: Definiert Leistungsgrenzen speicherbeschränkter Algorithmen auf GPUs neu
  2. Praktischer Wert: Bietet kritische technische Unterstützung für großskalige wissenschaftliche Berechnungen und KI-Anwendungen
  3. Ökosystem-Beitrag: Open-Source-Implementierung fördert Cross-Platform-GPU-Linear-Algebra-Ökosystem-Entwicklung

Anwendungsszenarien

  • Großskalige Matrix-SVD-Berechnungen
  • Wissenschaftliche Berechnungen und Ingenieur-Simulation
  • Dimensionsreduktion und Merkmalsextraktion im maschinellen Lernen
  • GPU-Anwendungen, die hochleistungsfähige lineare Algebra benötigen

Referenzen

Das Papier zitiert 86 verwandte Literaturquellen, die wichtige Arbeiten in GPU-Computing, linearen Algebra-Algorithmen, Julia-Sprachökosystem und anderen Bereichen abdecken und eine solide theoretische Grundlage für die Forschung bieten.