Mixed precision solvers with half-precision floating point numbers for Lattice QCD on A64FX processor

Dit onderzoek toont aan dat het gebruik van half-precisie (FP16) in gemengde precisie-oplossers voor rooster-QCD op de A64FX-processor, dankzij geïntroduceerde herschaalstappen, numeriek stabiel is en slechts een beperkte toename in iteraties vereist ten opzichte van dubbele precisie.

De kern

De Kern: Sneller Rekenen met een "Korte" Schaal

Stel je voor dat je een gigantisch, ingewikkeld puzzelraadsel moet oplossen. Dit raadsel is een simulatie van de kleinste bouwstenen van het universum (deeltjes zoals quarks en gluonen). Om dit op te lossen, gebruiken supercomputers wiskundige formules.

Normaal gesproken gebruiken deze computers een heel nauwkeurige manier van rekenen (zoals een liniaal met millimeters: FP64). Dit is precies, maar het kost veel tijd en energie om die millimeters te meten.

De onderzoekers van dit papier wilden een snellere manier vinden. Ze dachten: "Wat als we een liniaal gebruiken met alleen centimeters (FP16 of 'half-precisie')?"

• Het voordeel: Je kunt veel sneller meten omdat de getallen korter zijn. De computer kan meer berekeningen per seconde doen.
• Het probleem: Een centimeter-liniaal is niet nauwkeurig genoeg voor de allerlaatste stap van het raadsel. Als je te ver gaat, raken de getallen zelfs uit beeld (ze worden zo klein dat ze "verdwijnen" of underflow heten). Het raadsel wordt dan onoplosbaar of onstabiel.

De Oplossing: De "Verrekijker"-Truc

De onderzoekers (van RIKEN en andere Japanse instituten) hebben een slimme truc bedacht om dit op te lossen. Ze noemen het een gemengde precisie-oplosser.

Stel je voor dat je een lange reis maakt:

1. De Snelle Rit (FP16): Je rijdt eerst razendsnel over een snelweg met je "centimeter-liniaal". Je komt er heel snel, maar je weet niet precies waar je uitstapt.
2. De Controle (FP64): Aan het einde check je je positie met je nauwkeurige "millimeter-liniaal". Als je een beetje naast het doel staat, corrigeer je je koers.
3. De Herhaling: Je herhaalt dit proces.

Het probleem was: Bij het gebruik van de "centimeter-liniaal" (FP16) werden de getallen tijdens de rit zo klein dat ze verdwenen. Het was alsof je probeerde een druppel water te meten met een emmer; de druppel was te klein om te zien.

De Innovatie: Rescaling (Opnieuw Schalen)
De onderzoekers hebben een nieuwe methode bedacht die ze "rescaling" noemen.

• De Analogie: Stel je voor dat je een heel klein geluidje (een druppel water) moet horen in een luidrum. Als je het geluid te klein maakt, hoor je niets. Maar als je het geluid versterkt (opschalen) zodat het net binnen het bereik van je oren valt, kun je het horen. Zodra je het hebt verwerkt, zet je het volume weer terug.
• In de wiskunde betekent dit: Ze vermenigvuldigen de getallen tijdelijk met een groot getal zodat ze niet "verdwijnen" in de computer, rekenen ze af, en delen ze daarna weer door.

Wat hebben ze ontdekt?

1. Het werkt! Met hun nieuwe "versterk-truc" (rescaling) konden ze de snelle "centimeter-liniaal" (FP16) veilig gebruiken zonder dat de berekening instortte.
2. Het is snel: De berekeningen waren ongeveer drie keer zo snel als de oude, nauwkeurige methode en twee keer zo snel als de methode met "centimeters" die al bestond (FP32).
3. Minder fouten: Ze merkten op dat zonder hun truc, de computer veel "nuls" begon te zien (verdwijnende getallen). Met de truc bleven de getallen zichtbaar en nauwkeurig genoeg.

Waarom is dit belangrijk?

Deze studie is gedaan op de Fugaku, een van de krachtigste supercomputers ter wereld (in Japan).

• Voor de toekomst: De onderzoekers bereiden zich voor op de volgende generatie supercomputers (Fugaku NEXT), die nog meer krachtige "centimeter-rekenaars" (zoals GPU's met Tensor Cores) zullen hebben.
• Efficiëntie: Door deze truc te gebruiken, kunnen wetenschappers veel grotere en complexere simulaties van het universum draaien in minder tijd en met minder energie.

Samenvatting in één zin

De onderzoekers hebben een slimme manier gevonden om de "snelle maar onnauwkeurige" rekenkracht van moderne computers te gebruiken voor complexe natuurkunde, door de getallen tijdelijk op te rekken zodat ze niet verdwijnen, waardoor simulaties veel sneller gaan.

Technische samenvatting

Probleemstelling

Wetenschappelijke simulaties, zoals die in de Lattice QCD (Quantum Chromodynamica), vereisen doorgaans hoge rekennauwkeurigheid (FP64, double precision) om numerieke stabiliteit te garanderen. Echter, moderne hardware, zoals de A64FX-processor (gebruikt in de supercomputer Fugaku) en GPU's, biedt aanzienlijke prestatievoordelen bij het gebruik van lagere precisies, zoals FP32 (single) en FP16 (half-precision).

Hoewel mixed-precision solvers die FP32 combineren met FP64 al goed gevestigd zijn, is de toepassing van FP16 in deze context problematisch. De voornaamste uitdaging is de beperkte dynamische bereik van FP16 (exponentbereik van ongeveer $10^{-5}$ tot $6.5 \times 10^4$). Bij iteratieve methoden, zoals BiCGStab, neemt de norm van het residuvector tijdens de iteraties af. In FP16 leidt dit snel tot underflow (waarden die kleiner zijn dan het kleinste representeerbare getal en tot nul worden afgerond), wat de convergentie van de solver destabiliseert of volledig doet mislukken. Bestaande algoritmen die werken met FP32 presteren niet efficiënt of stabiel als ze direct worden overgezet naar FP16.

Methodologie

De auteurs onderzoeken de haalbaarheid van FP16 in mixed-precision lineaire solvers voor de Wilson-fermionmatrix (een vereenvoudigde formulering van de kwark-actie in Lattice QCD) op de A64FX-processor. De kern van hun aanpak is het introduceren van rescaling (opnieuw schalen) stappen om numerieke instabiliteit en underflow te voorkomen.

1. Algoritme-aanpassingen:

   • Iteratieve verfijning (Outer Loop): In de Richardson-iteratie (Algorithm 3) wordt het residuvector $r$ geschaald met een factor $\alpha$ zodat de norm een vooraf bepaald doelwaarde $s$ bereikt voordat het naar FP16 wordt omgezet.
   • BiCGStab Solver (Inner Loop): De auteurs passen een gerescaleerde versie van de BiCGStab-algoritme toe (Algorithm 4 en 5). Hierbij wordt het residuvector en de zoekrichtingsvectoren ($p, v$) continu geschaald met een factor $\gamma$ om te voorkomen dat waarden onder de ondergrens van FP16 zakken.
   • Oplossingsvector: In variatie Algorithm 5 wordt ook de oplossingsvector $x$ geschaald met een factor $\lambda$ om overflow te voorkomen als de oplossing veel groter is dan de invoer (bijvoorbeeld bij kleine eigenwaarden).
   • Herschaalfactoren: De schalingsfactoren ($\gamma$ en $\lambda$) worden dynamisch herberekend in FP32 (of FP64 voor globale reducties) om de nauwkeurigheid van de schaling te behouden, terwijl de zware rekenwerkzaamheden in FP16 plaatsvinden.

2. Implementatie op A64FX:

   • De implementatie maakt gebruik van de Scalable Vector Extension (SVE) van ARM, die 512-bit vectoren ondersteunt.
   • Er wordt gebruikgemaakt van het _Float16 datatype (ISO/IEC TC 18661-3:2015) in plaats van __fp16, omdat dit daadwerkelijke FP16-arithmetic ondersteunt in plaats van alleen opslag.
   • Voor sommaties en reducties wordt een hiërarchische aanpak gebruikt: sommatie binnen een thread in FP32, en globale reductie over threads en MPI-processen in FP64 om afrondingsfouten te minimaliseren.
   • De code is gebaseerd op de Bridge++ lattice QCD code-set, met specifieke optimalisaties voor de A64FX-architectuur (QXS-branch).

Belangrijkste Bijdragen

• Stabilisatie van FP16: Het paper introduceert een robuuste methode om FP16 te gebruiken in mixed-precision solvers voor Lattice QCD door systematische rescaling toe te passen in zowel de buitenste iteratieve verfijning als de binnenste BiCGStab-solver. Dit lost het underflow-probleem op dat eerder het gebruik van FP16 onmogelijk maakte.
• Prestatieverbetering: Het aantonen dat een naieve overdracht van FP32-algoritmen niet werkt, maar dat met de voorgestelde aanpassingen FP16 een praktische en snellere alternatief biedt.
• Hardware-geoptimaliseerde implementatie: Een efficiënte implementatie van FP16 op de A64FX-processor, inclusief het omzeilen van beperkingen in SIMD-instructies (zoals het ontbreken van compact-instructies voor 16-bit data) via expliciete conversies en opslagstrategieën.

Resultaten

De experimenten zijn uitgevoerd op de supercomputer Fugaku (RIKEN) met een rooster van $32^3 \times 64$ en een hopping parameter $\kappa = 0.13$.

• Convergentie:
  • Een standaard BiCGStab-solver in FP16 (zonder rescaling) faalde of convergeerde extreem traag (5518 matrix-vector vermenigvuldigingen) door instabiliteit.
  • Met de voorgestelde rescaling-methode (Algoritme 3 + 4 of 5) convergeerde de solver stabiel. Het aantal iteraties bleef binnen 20% van het aantal iteraties van de FP64-versie (bijv. ~850-920 iteraties versus 775 voor FP64).
• Snelheid en Prestaties:
  • Rekentijd: De totale doorlooptijd voor FP16 was ongeveer de helft van die van FP32 en ongeveer een derde van die van FP64.
  • FLOPS: De gemeten prestaties voor matrix-vector vermenigvuldiging waren:
    • FP64: ~2045 GFlops
    • FP32: ~3895 GFlops
    • FP16: ~8249 GFlops (een factor 4 hoger dan FP64, wat overeenkomt met de theoretische verwachtingen voor SVE).
  • Underflow: De analyse van het aantal nulwaarden in de vector toonde aan dat met de juiste rescaling (bijv. $s=128, \sigma=64$) underflow effectief werd onderdrukt, terwijl bij lagere schalingsfactoren onderflow optreedt en de convergentie vertraagt.

Betekenis en Toekomstperspectief

Dit werk is significant omdat het de weg vrijmaakt voor het gebruik van half-precisie (FP16) in complexe wetenschappelijke simulaties op toekomstige exascale-systemen, zoals de geplande "Fugaku NEXT" met NVIDIA Tensor Cores.

• Efficiëntie: Het demonstreert dat wetenschappelijke toepassingen die traditioneel als "FP64-only" werden beschouwd, kunnen profiteren van de enorme bandbreedte en rekenkracht van FP16 hardware zonder in te boeten aan nauwkeurigheid van het eindresultaat.
• Generaliseerbaarheid: De voorgestelde rescaling-methode is niet beperkt tot BiCGStab of Lattice QCD; het kan worden toegepast op andere mixed-precision preconditioners en iteratieve solvers om underflow te voorkomen.
• Toekomstig werk: De auteurs wijzen erop dat verdere studie nodig is voor complexere fermion-matrices (zoals Clover en Domain-wall) en voor het gebruik van BF16 (Brain Floating Point), evenals de implementatie op GPU's met Tensor Cores.

Samenvattend bewijst dit paper dat FP16 een praktisch en krachtig hulpmiddel is voor high-performance computing in de deeltjesfysica, mits de numerieke valkuilen van het beperkte bereik worden opgevangen door slimme algoritmische aanpassingen.