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 Algoritmo GPU-Residente Consapevole della Memoria per Accelerare la Bidiagonalizzazione di Matrici a Banda

Informazioni Fondamentali

  • ID Articolo: 2510.12705
  • Titolo: A GPU-resident Memory-Aware Algorithm for Accelerating Bidiagonalization of Banded Matrices
  • Autori: Evelyne Ringoot, Rabab Alomairy, Alan Edelman (MIT Computer Science & AI Laboratory)
  • Classificazione: cs.DC (Distributed, Parallel, and Cluster Computing), cs.MS (Mathematical Software)
  • Data di Pubblicazione: 14 ottobre 2025 (preprint arXiv)
  • Link Articolo: https://arxiv.org/abs/2510.12705

Riassunto

Questo articolo propone il primo algoritmo GPU-residente consapevole della memoria per accelerare la riduzione di matrici a banda a forma bidiagonale, un passaggio cruciale nella decomposizione ai valori singolari (SVD). Sebbene l'algoritmo possieda un elevato grado di parallelismo, è stato precedentemente considerato inadatto al calcolo GPU a causa della sua natura vincolata dalla larghezza di banda della memoria. Questa situazione è cambiata con l'evoluzione dell'hardware GPU, in particolare con l'aumento della memoria L1 per ogni streaming multiprocessor/unità di calcolo. Gli autori hanno basato il loro lavoro su precedenti algoritmi cache-efficienti bulge-chasing paralleli multi-core per CPU, ottimizzandoli per il throughput GPU. L'algoritmo è stato implementato come una singola funzione indipendente dall'hardware e dalla precisione dei dati su GPU NVIDIA, AMD, Intel e Apple Metal, supportando calcoli a mezza precisione, singola precisione e doppia precisione. Gli esperimenti dimostrano che l'algoritmo GPU supera le librerie CPU multi-thread ad alte prestazioni PLASMA e SLATE a partire da matrici 1024×1024, con miglioramenti di prestazioni superiori a 100 volte su matrici 32k×32k.

Contesto di Ricerca e Motivazione

Definizione del Problema

La decomposizione ai valori singolari (SVD) è uno strumento numerico fondamentale nella computazione scientifica, nell'apprendimento automatico e nell'analisi dei dati, con applicazioni diffuse nell'analisi delle componenti principali, nell'indicizzazione semantica latente, nell'approssimazione di basso rango e nel completamento di matrici. L'SVD su hardware moderno su larga scala segue tipicamente un processo a tre fasi:

  1. Riduzione da matrice densa a matrice a banda
  2. Riduzione da matrice a banda a matrice bidiagonale (bulge-chasing)
  3. Riduzione da matrice bidiagonale a matrice diagonale

Motivazione della Ricerca

Sebbene le implementazioni GPU della prima e della terza fase siano state ampiamente studiate, la seconda fase rimane insufficientemente esplorata su GPU moderne. Dongarra et al. nel 2014 hanno osservato che "gli acceleratori mostrano prestazioni scadenti nel trattare compiti computazionali a grana fine vincolati dalla memoria (come il bulge chasing), limitando i potenziali benefici dell'implementazione GPU della seconda fase".

Opportunità Tecniche

I recenti progressi nell'architettura GPU, in particolare:

  • Espansione della dimensione della cache L1 per streaming multiprocessor
  • Larghezza di banda on-chip superiore
  • Modelli di accesso alla memoria più flessibili
  • Migliore riutilizzo della cache e gerarchia di memoria più ampia

Questi miglioramenti hanno significativamente alterato l'equilibrio memoria-calcolo, creando nuove opportunità per la riprogettazione di algoritmi vincolati dalla memoria.

Contributi Principali

  1. Primo algoritmo bulge-chasing GPU-residente consapevole della memoria: Propone il primo algoritmo nativo GPU per la riduzione di matrici a banda a forma bidiagonale, sfruttando appieno le proprietà di memoria superiori delle GPU ad alte prestazioni di ultima generazione, con miglioramenti di prestazioni di 10-100 volte rispetto alle implementazioni CPU ottimizzate.
  2. Strategia bulge-chasing cache-efficiente per matrici a banda larga: Attraverso una nuova strategia cache-efficiente con riduzione progressiva della larghezza di banda, l'algoritmo GPU è in grado di gestire larghezze di banda più grandi rispetto al passato, alterando l'equilibrio ottimale della larghezza di banda tra la prima e la seconda fase in SVD su larga scala.
  3. Implementazione open-source indipendente dall'hardware e dalla precisione dei dati: Una libreria open-source basata su astrazioni del linguaggio Julia con una singola definizione di funzione che supporta calcoli a singola precisione, mezza precisione e doppia precisione su GPU NVIDIA, AMD, Intel e Apple.

Dettagli del Metodo

Definizione del Compito

Data una matrice a banda A ∈ ℝⁿˣⁿ con larghezza di banda BW, l'obiettivo è ridurla a forma bidiagonale B attraverso trasformazioni ortogonali, tale che A = UBVᵀ, dove U e V sono matrici ortogonali.

Architettura dell'Algoritmo

Flusso dell'Algoritmo Principale

Algoritmo 1: Riduzione da Matrice a Banda a Forma Bidiagonale 
             Utilizzando Vettori di Householder
Input: Larghezza di banda BW, larghezza piastrella interna TW, dimensione matrice n
1: for riduzione larghezza di banda i = (BW-1)/TW → 1 do
2:   Larghezza di banda target TBW = 1 + i·TW  
3:   for parallelo: ogni riga R = 1→n do
4:     bulge di riga k = R
5:     for ogni bulge di riga: j = 0, j+=1 do
6:       if 3(R-1) < j and k ≤ n then
7:         Calcola vettore HH per riga k per eliminare TW elementi 
            e applica alle righe sottostanti
8:         Calcola vettore HH per bulge di colonna generato più a sinistra 
            e applica alle colonne a destra
9:         Aggiorna valore k
10:      end if
11:    end for
12:  end for
13: end for

Implementazione del Kernel GPU Consapevole della Memoria

Algoritmo 2: Kernel GPU Consapevole della Memoria, TPB thread per blocco
1: Memoria thread: Ai (TW+1)
2: Memoria blocco: X (TW+1)  
3: Tutti i thread nel blocco collaborano: X ← A[k,..]
4: Sincronizza thread
5: Tutti i thread nel blocco collaborano: HH(X)
6: Tutti i thread nel blocco collaborano: A[k,..]← X
7: Sincronizza thread
8: for l: 0→(CBW + TW)/TPB - 1 do
9:   Thread i: calcola riga r = k + l·CPB + i
10:  Ai ← A[r, ...]
11:  HH(X, Ai)  
12:  A[r, ...]← Ai
13: end for

Punti di Innovazione Tecnica

1. Strategia di Partizionamento della Larghezza di Banda

Divide la matrice in piastrelle di larghezza di banda, eseguendo la riduzione per blocchi di larghezza di banda consecutivi, piuttosto che ridurre l'intera larghezza di banda in una volta. Questa strategia di partizionamento realizza:

  • Migliore riutilizzo della cache
  • Migliore località di memoria
  • Riduzione del sovraccarico di sincronizzazione

2. Pianificazione dei Blocchi Controllata

Introduce un parametro Max blocks che limita il numero di blocchi di thread concorrenti per unità di esecuzione GPU. Quando il numero di blocchi bulge-chasing richiesti dall'algoritmo supera il limite, adotta lo svolgimento del ciclo a livello software, con un singolo blocco assegnato a più compiti ed eseguiti sequenzialmente nel medesimo lancio del kernel.

3. Strategia di Separazione a Tre Cicli

Per garantire la correttezza e abilitare il parallelismo, l'algoritmo applica una separazione a tre cicli tra le scansioni di righe consecutive. Cioè, dopo il completamento di tre bulge di riga, la riga successiva può iniziare la scansione senza sovrapposizione degli accessi ai dati.

Configurazione Sperimentale

Ambiente Hardware

  • CPU: Intel Xeon Platinum 8462Y+ 32 core 64 thread
  • GPU: Include NVIDIA A100/H100/RTX4060, AMD MI250X/MI300X, Intel PVC 1100, Apple M1

Benchmark di Confronto

  • PLASMA (v25.5.27): Parallel Linear Algebra Software Package
  • SLATE (v2024.10.29): Scalable Linear Algebra Library
  • CUBLAS geam: Kernel di addizione matriciale come riferimento di larghezza di banda della memoria

Metriche di Valutazione

  • Tempo di esecuzione e rapporto di accelerazione
  • Throughput di memoria (DRAM, L1, L2)
  • Precisione numerica (errore relativo)
  • Utilizzo delle risorse hardware

Configurazione dei Test

  • Dimensione matrice: 1024×1024 a 65k×65k
  • Intervallo di larghezza di banda: 32 a 512
  • Precisione dei dati: FP16, FP32, FP64

Risultati Sperimentali

Risultati Principali di Prestazione

Confronto Prestazioni GPU vs Librerie CPU

  • Larghezza di banda piccola (32): Implementazione GPU 100-1000 volte più veloce di PLASMA, 10-300 volte più veloce di SLATE
  • Larghezza di banda grande (512): Versione GPU 2-9 volte più veloce di SLATE, 0,9-6 volte più veloce di PLASMA
  • Effetto della dimensione della matrice: Le prestazioni GPU crescono più velocemente dell'implementazione CPU con l'aumento della dimensione della matrice

Prestazioni su Hardware Diverso

  • NVIDIA H100: Prestazioni di riferimento
  • AMD MI300X: Circa 1,5-2 volte riduzione di prestazioni
  • Intel PVC: Circa 20 volte riduzione di prestazioni (nonostante cache più grande)
  • Apple M1: Buone prestazioni tra le GPU consumer

Risultati dell'Ottimizzazione degli Iperparametri

Scoperte chiave:

  1. La larghezza della piastrella interna è il fattore di prestazione più importante
    • Valore ottimale FP32: 32 (corrispondente a 128 byte di linea cache)
    • Valore ottimale FP64: 16 (corrispondente a 128 byte di linea cache)
  2. L'importanza dei parametri varia con la larghezza di banda:
    • Larghezza di banda 32: numero massimo di blocchi più critico
    • Larghezza di banda 128: numero di thread per blocco più critico

Analisi della Precisione Numerica

  • FP64: Errore prossimo alla precisione della macchina
  • FP32: Crescita dell'errore prevedibile dipendente dalla dimensione
  • FP16: Mantiene precisione accettabile entro matrici 16k×16k

Impatto dell'Evoluzione Hardware

  • MI300X mostra miglioramenti significativi rispetto a MI250X (cache L1 raddoppiata, aggiunto livello L2.5 unificato)
  • H100 continua a superare A100 (cache L1 aumentata del 33%, L2 aumentata del 25%)

Lavori Correlati

Sviluppo dei Risolutori SVD

  • Implementazioni CPU: PLASMA, SBR, FLAME, ELPA, LAPACK e altri focalizzati su trasformazioni Householder a blocchi e algoritmi bulge-chasing
  • Implementazioni GPU: Principalmente concentrate sulla prima fase (densa a banda) o SVD completo a una fase, evitando accessi alla memoria irregolari

Confronto con Risolutori di Autovalori

A differenza dei risolutori SVD, i risolutori a due fasi per problemi di autovalori a banda simmetrici hanno ampi lavori su sistemi ibridi CPU-GPU, ma la ricerca sulla riduzione da banda non simmetrica a bidiagonale rimane indietro.

Algoritmo Bulge-chasing

Originariamente proposto come algoritmo parallelo, ma gli studi sul caso bidiagonale su CPU sono estremamente rari. Al contrario, la ricerca su bulge-chasing da banda simmetrica a tridiagonale e forma di Hessenberg superiore è più diffusa.

Conclusioni e Discussione

Conclusioni Principali

  1. Superamento della barriera della larghezza di banda della memoria: Dimostra che gli algoritmi di algebra lineare vincolati dalla memoria non solo sono fattibili su GPU moderne, ma possono superare significativamente le prestazioni CPU
  2. Importanza delle caratteristiche hardware: La latenza della cache L1/L2 è più critica della dimensione, sottolineando l'importanza dell'accesso rapido alla memoria a bassa latenza
  3. Principi di progettazione dell'algoritmo: Le prestazioni ottimali derivano dalla progettazione algoritmica principiata, parametrizzazione consapevole della cache e infrastruttura di compilazione portabile

Limitazioni

  1. Modello di occupazione GPU: L'algoritmo ha un numero di blocchi GPU significativamente inferiore alle unità di esecuzione disponibili nella fase iniziale, richiedendo matrici sufficientemente grandi per utilizzare pienamente le risorse GPU
  2. Spillover di Registri: Larghezze di piastrella interna elevate possono causare spillover di registri, influenzando le prestazioni
  3. Dipendenza dall'Hardware: Variazioni significative di prestazioni tra diverse architetture GPU richiedono ottimizzazione specifica per hardware

Direzioni Future

  1. Ottimizzazione dell'Architettura: Sfruttare nuove funzionalità come CUDA 13.0 con spillover di memoria condivisa in registri
  2. Estensione dell'Algoritmo: Esplorare strategie di elaborazione per larghezze di banda più grandi
  3. Pipeline SVD Completa: Costruire una pipeline SVD completamente GPU-residente

Valutazione Approfondita

Punti di Forza

  1. Contributo Pioneristico: Prima implementazione di un algoritmo GPU-residente da banda a bidiagonale, colmando un importante vuoto tecnico
  2. Alto Valore Pratico: Miglioramenti significativi di prestazioni (10-1000 volte) conferiscono importante valore applicativo
  3. Innovazione Tecnica: Progettazione intelligente del partizionamento della larghezza di banda e consapevolezza della memoria adatta alle caratteristiche dell'architettura GPU
  4. Forte Portabilità: Implementazione single-source basata su Julia supporta GPU multi-vendor e multi-precisione

Carenze

  1. Analisi Teorica Insufficiente: Mancanza di analisi della complessità algoritmica e prove di convergenza
  2. Ambito di Test Limitato: Test principalmente su matrici sintetiche, mancanza di validazione in scenari di applicazione reale
  3. Complessità dell'Ottimizzazione dei Parametri: Richiede complessa ottimizzazione degli iperparametri per diversi hardware

Impatto

  1. Significato Accademico: Ridefinisce i confini di prestazione degli algoritmi vincolati dalla memoria su GPU
  2. Valore Pratico: Fornisce supporto tecnologico critico per il calcolo matriciale su larga scala e applicazioni AI
  3. Contributo all'Ecosistema: L'implementazione open-source promuove lo sviluppo dell'ecosistema di algebra lineare GPU multi-piattaforma

Scenari Applicabili

  • Calcolo SVD di matrici su larga scala
  • Computazione scientifica e simulazione ingegneristica
  • Ridimensionamento e estrazione di caratteristiche nell'apprendimento automatico
  • Applicazioni GPU che richiedono algebra lineare ad alte prestazioni

Bibliografia

L'articolo cita 86 lavori correlati, coprendo molteplici campi inclusi il calcolo GPU, gli algoritmi di algebra lineare e l'ecosistema del linguaggio Julia, fornendo una solida base teorica per la ricerca.