9.4 Stabil numerika és érzékenység

9.4.1 Helyzet és cél

Egy rendelésbecslő modellben a normalizáló log-likelihood tagot kell kiszámítanunk. A bemenet az egyes minták logitjainak vektora, például

$$\mathbf{z} = [1000, 1001, 1002]^T,$$

ahol minden $z_i$ dimenziótlan valós szám. A megfigyelési egység egyetlen mintapont, az eredmény pedig a skalár

$$L(\mathbf{z}) = \log \sum_{i=1}^{3} \exp(z_i).$$

Előzetes becslés: mivel a legnagyobb tag $\exp(1002)$ dominál, $L$ a pontos stabil alak $1002+\log(1+e^{-1}+e^{-2})\approx1002{,}4076$; a $1002+\log 3$ csak durva felső korlát. A naiv np.log(np.exp(z).sum()) megoldás azonban inf értéket adhat, mert $\exp(1000)$ már túlcsordul a float64 tartományon. Célunk, hogy ugyanezt a matematikai képletet stabil alakban számítsuk ki, és a véges kimenetet ellenőrizzük.

9.4.2 Miért csúszik el a számítás?

A float64 számábrázolás véges tartományú. A kerekítési hiba fogalmát már a Határérték mint közelítési folyamat alfejezetben bevezettük; most a túlcsordulást, alulcsordulást és kivonási hibát különítjük el.

Overflow (túlcsordulás): ha egy művelet eredménye nagyobb a legnagyobb ábrázolható számnál, a float64 inf-ként jelöli. Például $\exp(710) \approx 2.23 \cdot 10^{308}$, de $\exp(711)$ már túlcsordul. A math.exp(710) tipikus CPython-környezetben már OverflowError-t vált ki (míg a NumPy np.exp(710) jellemzően inf-et ad); a pontos viselkedés a művelettől és a könyvtártól függ; a pontos határ implementációfüggő, ezért ne a 710-re építsünk, hanem a stabil átírásra.

Underflow (alulcsordulás): ha az eredmény abszolút értéke kisebb a legkisebb ábrázolható pozitív számnál, az érték $0$-vá válik. Például $\exp(-750) < 10^{-324}$, ezért math.exp(-750) a gyakorlatban $0.0$. Ez önmagában nem mindig káros, de ha egy későbbi osztás nevezőjébe kerül, vagy ha logaritmust kellene rá alkalmazni, akkor hibát okoz.

Cancellation (kivonási hiba): két egymáshoz közeli, hosszú szám kivonásakor a közös vezető jegyek kiejtik egymást, és a megmaradó különbségben a bemenetek kerekítési hibája aránytalanul felnagyítódik. Klasszikus példa:

$$f(x) = \sqrt{x^2+1} - x.$$

Nagy $x$-re a két tag majdnem egyenlő; a kivonás elveszti az értékes jegyeket. A racionalizált alak viszont stabil:

$$f(x) = \frac{1}{\sqrt{x^2+1} + x}.$$

A két kifejezés algebrailag azonos, numerikusan azonban nagyon eltérő lehet. Ez a későbbi gyakori tévedés magja.

9.4.3 Relatív hiba és függvényérzékenység

A hiba szemléltetéséhez használjuk a relatív hibát. Ha $a$ a közelítő érték és $a_{\text{ref}}$ a referenciaérték, akkor

$$e_{\text{rel}} = \frac{|a - a_{\text{ref}}|}{|a_{\text{ref}}|}, \qquad a_{\text{ref}} \ne 0.$$

A relatív hiba dimenziótlan; például ha $a = 1002.4076$ és $a_{\text{ref}} = 1002.4076056$, akkor $e_{\text{rel}} \approx 5.6 \cdot 10^{-9}$.

A függvényérzékenység azt mutatja, hogy a bemeneti apró kerekítési hiba mennyire erősödik fel a kimenetben. A $\sqrt{x^2+1} - x$ esetén nagy $x$-nél a függvény kimenete kicsi, de a két bemeneti tag hibája nagy; a relatív hiba ezért robban. A stabil alakban már nincs ilyen kivonás, így a relatív hiba a float64 gépi pontosság közelében marad.

9.4.4 logsumexp: a stabil normalizáló logaritmus

A labor központi eszköze a log-sum-exp rövidítéssel ismert számítás. Az általános képlet $\mathbf{z} \in \mathbb{R}^n$ vektorra:

$$L(\mathbf{z}) = \log \sum_{i=1}^{n} \exp(z_i).$$

Ha a $z_i$ értékek nagyok, az exponenciálisak túlcsordulnak. A stabil alak kiemeli a legnagyobb értéket. Legyen

$$m = \max_i z_i,$$

ekkor minden $i$-re $z_i - m \le 0$, ezért $\exp(z_i - m) \le 1$, és nem lép fel túlcsordulás. A stabil képlet:

$$L(\mathbf{z}) = m + \log \sum_{i=1}^{n} \exp(z_i - m).$$

A helyesség levezetése:

$$ L = \log\Bigl(e^m \sum_i e^{z_i-m}\Bigr) = \log e^m + \log \sum_i e^{z_i-m} = m + \log \sum_i e^{z_i-m}. $$

A második tag legfeljebb $\log n$, tehát az eredmény mindig $m$ és $m+\log n$ közé esik.

Kézi kis kontroll. Legyen $\mathbf{z} = [10, 11, 12]$. Ekkor $m=12$, az eltérések $[-2, -1, 0]$, így

$$\sum_{i} e^{z_i - m} = e^{-2} + e^{-1} + e^{0} \approx 0.1353353 + 0.3678794 + 1 = 1.5032147.$$

A logaritmus $\log(1.5032147) \approx 0.4076056$, tehát

$$L = 12 + 0.4076056 \approx 12.4076056.$$

Ezt az eredményt használjuk ellenőrző invariantként: a stabil implementációnak ugyanezt kell adnia $\approx 10^{-12}$ tűréssel.

9.4.5 Skálázás mint előfeldolgozás

A logsumexp eltolás-invariáns: ha minden elemhez ugyanazt a $c$ konstanst adjuk, az eredmény $c$-vel eltolódik, mert

$$L(\mathbf{z} + c\mathbf{1}) = L(\mathbf{z}) + c.$$

Ez az alapja a max-központosításnak. A skálázás ezzel szemben nem invariáns: ha minden elemet ugyanazzal a $s$ skalárral szorzunk, a $\exp(s z_i)$ szerkezet miatt az eredmény általában megváltozik, és nem egyszerű eltolással tér el. A softmax például a skálázásra hőmérséklet-értelemben változik:

$$p_i = \frac{\exp(z_i)}{\sum_j \exp(z_j)}, \qquad p_i^{(\tau)} = \frac{\exp(z_i/\tau)}{\sum_j \exp(z_j/\tau)},$$

ahol $\tau > 0$ dimenziótlan hőmérséklet. A $\tau$ növelése az eloszlást laposabbá, csökkentése csúcsosabbá teszi. Numerikai szempontból a skálázás segíthet az overflow elkerülésében, de megváltoztatja a modell szemantikáját; ezért a kísérleti jegyzőkönyvben rögzíteni kell, hogy milyen skálán dolgoztunk.

9.4.6 Teljes munkapélda: naiv és stabil log-összeg-exp

Az alábbi kód a Northstar kosárlogitokhoz hasonló méretű nagy bemeneten hasonlítja össze a két alakot. A np.logaddexp a NumPy beépített stabil logsumexp-alapművelete.

import math
import numpy as np

# Bemenet: egy mintapont logitvektora, shape (3,), dtype float64.
z_large = np.array([1000.0, 1001.0, 1002.0])

def naive_logsumexp(z):
    # Mathematikai képlet, numerikusan törékeny.
    return np.log(np.sum(np.exp(z)))

def stable_logsumexp(z):
    # Stabil képlet: kiemeljük a maximumot.
    m = np.max(z)                      # skalár
    return m + np.log(np.sum(np.exp(z - m)))

for name, fn in [("naiv", naive_logsumexp), ("stabil", stable_logsumexp)]:
    value = fn(z_large)
    print(f"{name:6s} L = {value!r}, véges: {np.isfinite(value)}")

# Kézi kontroll a kis [10, 11, 12] bemeneten.
z_small = np.array([10.0, 11.0, 12.0])
expected = 12.0 + math.log(math.exp(-2.0) + math.exp(-1.0) + 1.0)
got = stable_logsumexp(z_small)
print(f"kézi várt: {expected:.12f}")
print(f"stabil   : {got:.12f}")
print(f"hiba     : {abs(got - expected):.3e}")

# NumPy beépített változat ugyanarra a kis bemenetre.
np_check = np.logaddexp.reduce(z_small)
print(f"np.logaddexp.reduce: {np_check:.12f}")

A nem futtatott kimenet szemléltető; a várható értékek elméleti alapon:

naiv   L = inf, véges: False
stabil L = 1002.4076055676442, véges: True
kézi várt: 12.407605567644
stabil   : 12.407605567644
hiba     : ~1e-15
np.logaddexp.reduce: 12.407605567644

A naiv kimenet a float64 overflow miatt inf, ezért a végesérték-ellenőrzés (np.isfinite) hamis lesz. A stabil változat véges, és a kis bemeneten a kézi kontrollal egyezik.

9.4.7 Eredmény értelmezése és ellenőrző invariantok

A stabil L értéke mindig $[m, m+\log n]$ intervallumba esik. A fenti nagy bemeneten $m = 1002$ és $\log 3 \approx 1.0986$, így $L \in [1002, 1003.0986]$; a kapott $1002.4076$ ezen belül van. Ez jó önellenőrzés: ha az eredmény kívül esik az intervallumon, akkor vagy a max-számítás hibás, vagy az összegzésben csúszott be inf/nan.

A NumPy np.logaddexp.reduce(z_small) ugyanazt a stabil számítást valósítja meg; a kis bemeneten a két implementáció közötti különbség a float64 kerekítéséből adódik, és $10^{-13}$ alatt marad. Tömbalakok: bemenet (n,), kimenet skalár; kétdimenziós bemenetnél axis=1 esetén a kimenet (rows,) lesz, mert az összegzés a tengely mentén történik.

9.4.8 Feladatok

1. Értelmezés kézi számolással. Számítsd ki kézzel a stabil logsumexp értékét a $\mathbf{z} = [-2, -1, 0]$ bemenetre. Add meg $m$ értékét, az eltéréseket, az összeget, a logaritmust és a végső $L$-t. Miért nincs itt sem overflow, sem underflow?

2. Implementálás saját függvénnyel. Írj stable_logsumexp(z, axis=None) függvényt, amely egy- és kétdimenziós np.ndarray bemeneten is működik. A kétdimenziós esetben axis=1 mentén számolj, és a kimenet shape-je legyen (rows,). Ellenőrizd a

   Z = np.array([[0.0, 1.0], [1000.0, 1001.0]])

   mátrixszal; a várt kimenet elméletileg [1.31326169, 1001.31326169] körül alakul. Hasonlítsd össze a np.logaddexp.reduce(Z, axis=1) eredményével, és írasd ki a maximális abszolút eltérést.

3. Hibakeresés cancellation-nel. Az alábbi függvény nagy x-re hibásan 0.0-t adhat, pedig matematikailag pozitív értéket várunk:

   def log_diff(x):
       return np.log(x + 1.0) - np.log(x)

   1. Mi a hiba neve? (b) Írd át a kifejezést stabil alakra a log1p használatával. (c) Ellenőrizd x = 1e16 esetén: a naiv 0.0, a stabil pedig körülbelül 1e-16; magyarázd, miért.

4. Megváltozott feltétel: cent és euró. Egy tároló a kosárértéket EUR-centben adja meg: basket_cents = np.array([99999995.0, 100000005.0, 100000015.0]). (a) A naiv softmax np.exp(z) / np.exp(z).sum() miért válik nan-ná? (b) Skálázd át a bemenetet euróra z_eur = basket_cents / 100.0, és számítsd ki a softmaxot; ezután számítsd ki a cent-alapú stabil softmaxot z_cent - z_cent.max() értékekkel. Hasonlítsd össze a két eredményt: mit mutat ez a skálázás szemantikai hatásáról? A cél az eredeti cent-logitok valószínűség-arányainak megőrzése; melyik megoldás teszi ezt?

9.4.9 Megoldások és indoklás

1. Megoldás. $\mathbf{z} = [-2, -1, 0]$, tehát $m = \max_i z_i = 0$. Az eltérések $[-2, -1, 0]$. Az exponenciálisok: $e^{-2} \approx 0.1353353$, $e^{-1} \approx 0.3678794$, $e^{0} = 1$. Összegük $1.5032147$. A logaritmus $\log(1.5032147) \approx 0.4076056$, így $L = 0 + 0.4076056 \approx 0.4076056$. Azért nincs sem overflow, sem underflow, mert minden eltérés $-2$ és $0$ közé esik, az exponenciálisak $0.135$ és $1$ között maradnak, amelyek a float64 normál tartomány kényelmes belsejében vannak.

2. Megoldás. Egy lehetséges implementáció:

   def stable_logsumexp(z, axis=None):
       z = np.asarray(z, dtype=np.float64)
       m = np.max(z, axis=axis)
       if axis is None:
           return m + np.log(np.sum(np.exp(z - m)))
       m_expanded = np.expand_dims(m, axis=axis)
       return m + np.log(np.sum(np.exp(z - m_expanded), axis=axis))
   
   Z = np.array([[0.0, 1.0], [1000.0, 1001.0]])
   print(stable_logsumexp(Z, axis=1))
   print(np.logaddexp.reduce(Z, axis=1))
   print(np.max(np.abs(stable_logsumexp(Z, axis=1) -
                       np.logaddexp.reduce(Z, axis=1))))

   A várt kimenet elméletileg:

   [   1.31326169 1001.31326169]
   [   1.31326169 1001.31326169]
   ~2e-15

   Az axis=1 esetén a kimenet shape (2,), mert a soronkénti maximum plusz a soronkénti log-összeg vektoros művelet. A np.logaddexp.reduce ugyanazt a stabil számítást végzi; az eltérés a gépi kerekítésből adódik, és $10^{-14}$ nagyságrendű marad.

3. Megoldás. (a) A hiba neve cancellation vagy kivonási hiba: $\log(x+1)$ és $\log(x)$ nagy $x$-re egymáshoz nagyon közeli számok, a float64 pedig a kivonáskor elveszti a hasznos jegyeket. (b) A stabil átírás:

   def log_diff_stabil(x):
       return np.log1p(1.0 / x)

   mert $\log(x+1)-\log(x) = \log(1 + 1/x)$. (c) $x = 10^{16}$ esetén a stabil alak log1p(1e-16) értéke körülbelül $10^{-16}$. A naiv log(1e16 + 1) - log(1e16) azért adhat 0.0-t, mert a float64 kerekítésében $10^{16}+1$ pontosan $10^{16}$-ként tárolódik, így a kivonás valódi különbsége elvész, mielőtt a logaritmus megtörténhetne. A stabil alakban az $1/x$ hányados $10^{-16}$-ként megmarad, és a log1p a kis argumentumot pontosan kezeli.

4. Megoldás. (a) A centek $\exp(10^8)$ értéket adnak, ami túlcsordul a float64-ben: inf. Az inf/inf hányados nan, ezért a softmax nan lesz. (b) Az euró-alapú softmax:

   z_eur = basket_cents / 100.0
   p_eur = np.exp(z_eur) / np.exp(z_eur).sum()
   # elméletileg: [0.300606, 0.332229, 0.367165]

   A cent-alapú stabil softmax:

   z_st = basket_cents - basket_cents.max()
   p_cent = np.exp(z_st) / np.exp(z_st).sum()
   # elméletileg: [2.06e-9, 4.54e-5, 9.999546e-1]

   A két eredmény eltérő, mert a skálázás nem pusztán numerikus trükk, hanem a hőmérsékletet módosítja: a cent-skáláról euróra osztás a logitokat $100$-zal osztja, így a valószínűségek laposabbá válnak. Ha az eredeti cent-logitok arányait akarjuk megtartani, a max-eltolásos cent-alapú stabilitás a helyes választás; az euróra osztás más modellt ad, amit a kísérleti jegyzőkönyvben explicit jelezni kell.

9.4.10 Gyakori tévedés és önellenőrzés

A leggyakoribb tévedés, hogy algebrailag azonos képleteket numerikusan is felcserélhetőnek tekintünk. Ellenpélda a sqrt(x**2+1)-x és az $1/(\sqrt{x^2+1}+x)$ pár: nagy $x$-nél az első akár teljesen rossz, a második pontos marad. Itt kapcsolódunk vissza a Határérték mint közelítési folyamat kerekítésihiba-tanulságához: a véges pontosság miatt a lépésköz csökkentése vagy a naiv képlet egyszerűsége nem mindig javít. A Vektorizálás és mért teljesítmény mérési fegyelme pedig arra emlékeztet, hogy a numerikus egyezést ne szemre, hanem np.allclose-zal és rögzített tűréssel ellenőrizzük. Az II. ellenőrző pont: adatból ellenőrzött feature-tábla feature-táblában ezt az elvet követtük: kontrollösszegek és toleranciás összehasonlítás.

Önellenőrző lista a laborhoz:

• Véges-e a kimenet? np.isfinite a nagy bemeneten is igaz.
• Illeszkedik-e az eredmény az $[m, m+\log n]$ intervallumba?
• A kis $[10,11,12]$ bemeneten a stabil érték $12.407605567644$ közelében van-e?
• A np.logaddexp.reduce eredménye $10^{-13}$ tűrésen belül egyezik-e a kézi függvénnyel?
• A cent-alapú és euró-alapú softmax különbségét megértettük-e, és a jegyzőkönyvben rögzítettük-e a skálaválasztást?

A logsumexp pontos értékét a maximum kiemelésével kell kiszámítani; a max(z) + log(n) csak felső korlát. A skalár math.exp túlcsordulási kivétele és a NumPy tömbműveletében keletkező inf nem azonos hibajelenség, ezért eltérő ellenőrzést igényelhetnek.

9.4.11 Rövid összefoglaló

A stabil numerika azt jelenti, hogy ugyanazt a matematikai célt olyan alakban számoljuk, amely a float64 véges tartományán belül marad. Az overflow inf-et, az underflow hamis 0-t, a cancellation pedig értékes jegyek elvesztését okozza. A logsumexp stabil változata a maximum kiemelésével kerüli el az overflow-t; a kézi $[10,11,12]$ kontroll $12.407605567644$-hez közeli eredményt ad. A skálázás nem ártalmatlan: a softmax hőmérsékletét, így a modell jelentését is megváltoztatja. A következő alkalmazás a Tanítás, validáció és érintetlen teszt prediktív pipeline lesz, ahol a log-loss stabil számítása és a softmax kalibrációja ugyanezt a fegyelmet igényli.