Sneller floating-point rekenen met de nieuwe API van Rust

Floating-point berekeningen zijn vaak trager dan integer-berekeningen omdat de compiler conservatief is in hoe hij je code optimaliseert. Hoewel sommige programmeertalen al een vorm van oplossing hadden, had Rust tot nu toe geen goede, stabiele manier om met deze beperking om te gaan. Vanaf versie 1.98 staat Rust het nu echter toe om de compiler te vertellen dat hij de code verder mag optimaliseren—maar mét extra controle, zodat je nog steeds numerieke algoritmen kunt schrijven met minimale afrondingsfouten.

In dit artikel leer je:

  • Waarom de compiler standaard floating-point berekeningen minder optimaliseert dan integer-berekeningen.
  • De nieuwe API van Rust om deze beperking op te lossen.
  • Voorbeelden van het gebruik van deze nieuwe API, de impact op de snelheid en hoe je kunt beheren waar dit wordt toegepast.

Het optellen van integers is snel

Ik begin met een voorbeeld met integers als referentiepunt voor wat er op het gebied van prestaties mogelijk is. Om de snelste codegeneratie te krijgen, vertel ik Rust dat we niet meer in 2004 leven en dat hij CPU-instructies mag genereren die moderne hardware vereisen (namelijk x86-64 machines van de afgelopen 10 jaar). Specifiek wordt alle code in dit artikel gecompileerd met RUSTFLAGS="-C target-cpu=x86-64-v3". (Voor maximale compatibiliteit zou je in een echte implementatie een fallback-versie voor oudere computers kunnen bieden.)

Hier is een Rust-functie om een slice van int64-getallen op te tellen:

fn naive_sum_i64(values: &[i64]) -> i64 {
    let mut total = 0;
    for value in values {
        total += value;
    }
    total
}

Om dit te benchmarken, maak ik een array van integers aan in NumPy:

import numpy as np
DATA_INT = np.ones((1_000_000,), dtype=np.int64)
assert naive_sum_i64(DATA_INT) == 1_000_000

En nu kan ik de snelheid van het optellen van deze array meten:

CodeVerstreken µ-secondenCPU-instructies per waarde
naivesumi64(DATA_INT)168.10.5

Dat is 0,5 CPU-instructie per waarde! Hoe kan dat? Waarschijnlijk gebruikt de compiler gespecialiseerde Single Instruction, Multiple Data (SIMD) CPU-instructies, die batch-operaties uitvoeren op meerdere waarden tegelijk. De i7-12700K CPU die ik hier gebruik heeft 256-bit SIMD-instructies, wat betekent dat hij specifieke operaties kan uitvoeren op vier 64-bit integers tegelijk. Als er een gespecialiseerde SIMD-sommatie instructie is, hoeft de CPU slechts 250.000 keer te loopen en in elke iteratie 4 integers op te tellen.

En dat is inderdaad het geval:

CodeVerstreken µ-secondenCPU-instructies256-bit SIMD integer instructies
naivesumi64(DATA_INT)156.8521,280250,003

Kortom, door een gespecialiseerde SIMD-instructie te gebruiken, kan mijn CPU integers zeer snel optellen.

Is het optellen van floats traag?

Maar hoe zit het met floats—zijn die ook zo snel? Opnieuw maak ik een miljoen floating-point waarden aan:

# Array van 1M float64 waarden tussen 0 en 1.
DATA = np.random.random((1_000_000,))

Ik implementeer een eenvoudige floating-point sommatie-functie:

fn naive_sum(values: &[f64]) -> f64 {
    let mut total = 0.0;
    for value in values {
        total += value;
    }
    total
}

En vergelijk de prestaties van het optellen van integers en floats:

CodeVerstreken µ-secondenCPU-instructies256-bit SIMD integer instructies256-bit SIMD float instructies
naivesumi64(DATA_INT)151.9521,214250,0030
naive_sum(DATA)595.21,458,26900

De floating-point sommatie is veel trager dan de integer-sommatie, en de compiler heeft geen SIMD-float operaties gebruikt. Waarom dit verschil?

Floating-point bewerkingen zijn niet associatief

Net als de meeste compilers optimaliseert Rust je code in release-modus door deze op verschillende manieren te transformeren om hem (hopelijk) sneller te maken. Maar compilers doen dit met een belofte: de geoptimaliseerde code zal zich precies hetzelfde gedragen als de ongeoptimaliseerde code.

Als ik drie integers $a, b,$ en $c$ optel, geldt dat $a + (b + c) == (a + b) + c$. Dit geeft de compiler veel ruimte om te optimaliseren, bijvoorbeeld door SIMD-operaties te gebruiken die de volgorde van de optellingen licht kunnen veranderen.

Floating-point getallen zijn anders. Omdat floating-point getallen een enorm bereik beslaan, van minuscuul tot gigantisch, resulteert het optellen van een zeer groot getal bij een zeer klein getal in datzelfde grote getal:

print(
    "Doet het toevoegen van een klein getal niets?",
    1e16 + 1.0 == 1e16
)
# Output: Doet het toevoegen van een klein getal niets? True

Breder gesteld is voor floating-point getallen $a + (b + c)$ niet altijd hetzelfde als $(a + b) + c$, althans zodra je meerdere getallen achter elkaar optelt. Stel dat ik een array heb die begint met 1e16 gevolgd door vele 1.0-waarden, en een andere die het omgekeerde is. Het optellen van deze arrays zal verschillende resultaten geven:

import math
HIGH_VALUE_FIRST = np.ones((1_000_000,), dtype=np.float64)
HIGH_VALUE_FIRST[0] = 1e16
HIGH_VALUE_LAST = np.ones((1_000_000,), dtype=np.float64)
HIGH_VALUE_LAST[-1] = 1e16

print(
    "Is de som hetzelfde?",
    naive_sum(HIGH_VALUE_FIRST) == naive_sum(HIGH_VALUE_LAST)
)
# Output: Is de som hetzelfde? False

Omdat de volgorde van optellen invloed heeft op het resultaat, gaat de compiler er (terecht) vanuit dat ik om een specifieke reden deze volgorde heb gevraagd. Als gevolg hiervan zal de compiler deze operaties niet herordenen, noch enige andere optimalisatie toepassen die de resultaten zou kunnen veranderen, zelfs als de resulterende code trager is.

De nieuwe algebraïsche operatoren van Rust: de compiler vertellen wanneer hij flexibel mag zijn

Hoewel conservatisme vanuit de compiler de juiste standaardinstelling is, weet jij als programmeur soms dat het herordenen van bewerkingen geen probleem is. In die situatie is het handig om de compiler te kunnen vertellen dat hij op de ene plek de volgorde niet mag veranderen, maar op een andere plek wel.

Vanaf Rust 1.98 is er een nieuwe functie die dit mogelijk maakt. Naast de normale rekenkundige bewerkingen voor floating-point getallen, is er een nieuwe set zogenaamde "algebraïsche" rekenkundige operatoren. Volgens de documentatie "staan deze de compiler toe om floating-point operaties te optimaliseren met gebruik van alle gebruikelijke algebraïsche eigenschappen van reële getallen", inclusief het wijzigen van de volgorde van bewerkingen.

De releasedatum van Rust 1.98 is 20 augustus 2026. Aangezien ik dit artikel daarvoor schreef, heb ik de code uitgevoerd met het "beta"-kanaal, dat dezelfde functionaliteit bevat.

Een voorbeeld: Geoptimaliseerde pairwise-sommatie

Laten we deze operatoren in actie zien. Omdat het optellen van floating-points onverwacht gedrag kan vertonen (denk aan 1e16 + 1.0 == 1e16), zijn er verschillende algoritmen die je kunt gebruiken om afrondingsfouten bij grote hoeveelheden getallen te minimaliseren. De afweging tussen deze algoritmen ligt meestal bij de snelheid versus de nauwkeurigheid.

numpy.sum() maakt voornamelijk gebruik van een algoritme genaamd pairwise summation, dat een goede balans biedt tussen snelheid en het verminderen van geaccumuleerde fouten. Het heeft zeker minder foutaccumulatie dan het simpelweg één voor één optellen in volgorde, zoals ik deed met naive_sum().

Het basisidee van pairwise-sommatie is om de array in tweeën te splitsen, elke zijde recursief op te tellen met hetzelfde algoritme, en vervolgens de twee resulterende floating-point getallen op te tellen. Onder een bepaalde drempelwaarde voor de array-grootte (128 in het geval van NumPy) gebeurt de sommatie normaal... en dat kan in elke willekeurige volgorde.

Laten we dit algoritme implementeren in Rust. Bij het optellen van de twee bovenste floating-point getallen gebruik ik normale optelling, zodat de compiler de operaties niet kan herordenen. Zodra de functie de drempel bereikt voor normale sommatie, schakel ik over op algebraïsche optellingen, omdat ik op dat punt niet om de volgorde geef, maar alleen om snelheid.

fn pairwise_sum(values: &[f64]) -> f64 {
    let n = values.len();
    if n > 128 {
        // Nauwkeurige optelling van twee recursieve toepassingen
        // van het algoritme:
        let half = n / 2;
        pairwise_sum(&values[0..half]) + pairwise_sum(&values[half..n])
    } else {
        // 😎 Normale optelling, waarbij de volgorde niet uitmaakt
        // voor het algoritme. Daarom is een algebraïsche optelling
        // prima—en dat geeft de compiler toestemming om
        // agressief te optimaliseren.
        let mut total: f64 = 0.0;
        for value in values {
            total = total.algebraic_add(*value);
        }
        total
    }
}

Deze implementatie van pairwise-sommatie is snel; hij is zelfs sneller dan de implementatie van NumPy:

CodeVerstreken µ-secondenCPU-instructies256-bit SIMD float instructies
naive_sum(DATA)563.11,458,2790
np.sum(DATA)190.72,191,7670
pairwise_sum(DATA)144.51,298,028270,336

En het heeft dezelfde (algemene) precisie als de implementatie van NumPy, aangezien het hetzelfde algoritme hanteert:

from math import fsum
# fsum() gebruikt een nauwkeurig sommatie-algoritme dat elke
# vermijdbare floating-point fout elimineert.
assert fsum(HIGH_VALUE_FIRST) == fsum(HIGH_VALUE_LAST)
correct_sum = fsum(HIGH_VALUE_FIRST)

print("naive_sum() error:   ", naive_sum(HIGH_VALUE_FIRST) - correct_sum)
print("np.sum() error:      ", np.sum(HIGH_VALUE_FIRST) - correct_sum)
print("pairwise_sum() error:", pairwise_sum(HIGH_VALUE_FIRST) - correct_sum)

# Output:
# naive_sum() error:    -1000000.0
# np.sum() error:       -14.0
# pairwise_sum() error: -6.0

Het verschil in fout tussen np.sum() en pairwisesum() is niet significant, dat is simpelweg geluk; het belangrijkste punt is dat ze een vergelijkbaar lage orde van grootte aan fouten hebben vergeleken met naivesum(). Als extra bonus is de codestructuur die pairwise_sum() gebruikt om SIMD-code te stimuleren veel beknopter dan die van NumPy.

Een ander voorbeeld: Som van gekwadrateerde verschillen (SSD)

Je kunt meer doen dan alleen getallen optellen met algebraïsche operatoren. In de volgende voorbeelden bereken ik de som van gekwadrateerde verschillen (sum of squared differences of SSD) tussen twee arrays.

Ik heb hiervoor twee arrays nodig:

DATA1 = np.random.random((1_000_000,))
DATA2 = np.random.random((1_000_000,))

Hier is een implementatie met normale rekenkundige operatoren:

fn ssd_normal(arr1: &[f64], arr2: &[f64]) -> f64 {
    assert_eq!(arr1.len(), arr2.len());
    let mut total = 0.0;
    for (val1, val2) in arr1.iter().zip(arr2) {
        total += (val1 - val2).powi(2);
    }
    total
}

En hier is hoe het eruitziet met algebraïsche operatoren; het is voor mij niet duidelijk of f64.powi() algebraïsche operatoren zal gebruiken, dus heb ik dit expliciet gedaan:

fn ssd_optimized(arr1: &[f64], arr2: &[f64]) -> f64 {
    assert_eq!(arr1.len(), arr2.len());
    let mut total: f64 = 0.0;
    for (val1, val2) in arr1.iter().zip(arr2) {
        // 😎 Algebraïsche operaties om meer compiler-optimalisaties toe te staan:
        let diff = val1.algebraic_sub(*val2);
        let squared_diff = diff.algebraic_mul(diff);
        total = total.algebraic_add(squared_diff);
    }
    total
}

Eerst een snelle test om te controleren of de resultaten vergelijkbaar zijn:

print(ssd_normal(DATA1, DATA2))      # 166770.0055951995
print(ssd_optimized(DATA1, DATA2))   # 166770.00559520238

Vervolgens meet ik de prestaties:

CodeVerstreken µ-secondenCPU-instructies per waarde
ssd_normal(DATA1, DATA2)628.74.5
ssd_optimized(DATA1, DATA2)371.11.0

Door algebraïsche operaties te gebruiken, kan de compiler code genereren die op mijn computer twee keer zo snel draait.

Versnel je numerieke code!

Het pairwise-sommatie algoritme is een uitstekend voorbeeld van waarom je beide soorten operatoren wilt hebben: strikte en lakse.

  • Als je alleen strikte, in-volgorde operaties gebruikt, is de code trager.
  • Als je alleen lakse algebraïsche "optimaliseer-zo-als-je-wilt" operaties gebruikt, zou de compiler het algoritme zelf weg kunnen optimaliseren, waardoor de nauwkeurigheid die het algoritme beoogt verloren gaat.

De implementatie die ik hierboven liet zien profiteert van het gebruik van zowel normale optelling (voor nauwkeurigheid) als algebraïsche optelling (voor snelheid) in verschillende delen van het algoritme.

Als je numerieke code schrijft met Rust, kan jouw code hier waarschijnlijk ook van profiteren—probeer het eens zodra je Rust 1.98 kunt gebruiken. En als je Rust nog niet gebruikt, is dit een reden meer om over te stappen.