Aritmetica di macchina - floating point, precisione di macchina e cancellazione

Prerequisiti: come si scrivono i bit di un numero in virgola mobile è in Numeri in virgola mobile IEEE 754Virgola fissa e virgola mobile; formato IEEE 754 a 32 e 64 bit (segno, esponente in eccesso 127, mantissa con 1 implicito); conversioni svolte nei due versi; valori speciali, denormalizzati, intervallo, precisione e arrotondamento.Numeri in virgola mobile IEEE 754 →; qui interessa invece quanto sbaglia il calcolatore quando fa un conto e come si progetta una formula che sbaglia poco. Il calcolo numerico nasce da questo fatto: il calcolatore non lavora con i numeri reali, ma con un insieme finito di numeri, e ogni risultato porta con sé un errore.

Struttura dell'elaboratore: quello che serve al calcolo numerico

Il dettaglio dell'hardware è in Architettura e organizzazione di un calcolatoreDifferenza tra architettura (ciò che vede il programmatore) e organizzazione (come è realizzata); struttura e funzione; le quattro funzioni e i quattro componenti; macchina di von Neumann e IAS; generazioni tecnologiche e legge di Moore.Architettura e organizzazione di un calcolatore →; per il calcolo numerico bastano tre fatti.

• Modello di von Neumann. La CPU (unità di controllo, ALU che esegue le operazioni, registri) legge dalla memoria centrale sia le istruzioni del programma sia i dati, li elabora e riscrive i risultati in memoria. Ogni dato occupa un numero fissato di bit: da qui nasce l'insieme finito di numeri di questa nota.
• Unità di informazione. Un bit vale $0$ o $1$; un byte sono $8$ bit; i processori attuali lavorano su parole di $64$ bit. Con $n$ bit si hanno $2^n$ configurazioni diverse, quindi al più $2^n$ numeri rappresentabili.
• Interi. Si usano in complemento a due (Rappresentazione dei numeri interi con segnoInteri con segno su n bit: modulo e segno, complemento a 1, complemento a 2 ed eccesso K; intervalli rappresentabili, calcolo dell'opposto, estensione del segno.Rappresentazione dei numeri interi con segno →): con $n$ bit l'intervallo è $[-2^{n-1},\,2^{n-1}-1]$. Le operazioni tra interi sono esatte finché il risultato resta nell'intervallo; fuori si ha overflow. Esempio a $8$ bit: l'intervallo è $[-128,\,127]$; $-5$ si scrive invertendo i bit di $5=00000101$ e sommando $1$: $11111010+1=11111011$ (controllo: $-128+64+32+16+8+2+1=-5$ ✓); $127+1$ dà $10000000=-128$ (il risultato «si avvolge»). In MATLAB gli interi int8, int32 invece saturano (int8(127)+1 vale 127); in Python gli interi non hanno limite.
• Reali. Si usa la virgola mobile in $64$ bit (doppia precisione: $1$ bit di segno, $11$ di esponente, $52$ di mantissa), che è l'argomento del resto della nota. È il tipo predefinito di MATLAB e il float di Python.
• Costo. Il tempo di un algoritmo numerico si misura in flop (operazioni in virgola mobile). Un processore da tavolo esegue circa $10^9$–$10^{11}$ flop/s, ma leggere un dato dalla memoria centrale costa molto più di un'operazione: per questo la memoria è a livelli (registri, cache, RAM, Gerarchia di memoria e principio di localitàCaratteristiche delle memorie (posizione, capacità, unità di trasferimento, metodo di accesso, prestazioni, volatilità); compromesso costo-capacità-velocità e gerarchia a livelli; località temporale e spaziale con esempi; hit, miss e tempo medio di accesso con calcolo svolto.Gerarchia di memoria e principio di località →). Esempio: risolvere un sistema lineare $1000\times1000$ richiede circa $\frac23\cdot10^9$ flop (Fattorizzazione LU, pivoting e CholeskyRisolvere $A\mathbf x=\mathbf b$ con Gauss equivale a fattorizzare $A=LU$ ($L$ triangolare inferiore con $1$ sulla diagonale e i moltiplicatori $l_{ik}=\frac{a_{ik}^{(k)}}{a_{kk}^{(k)}}$ sotto la diagonale, $U$ la matrice a scala finale) e poi risolvere $L\mathbf y=\mathbf b$ (in avanti) e $U\mathbf x=\mathbf y$ (all'indietro). Esiste (ed è unica) se i minori principali di testa $\det A_k$, $k<n$, sono $\ne0$. Costo $\frac23n^3$ flop per la fattorizzazione, $2n^2$ per ogni termine noto; $\det A=\prod u_{kk}$. Con un pivot piccolo l'algoritmo è instabile: il pivoting parziale (scambio di righe che porta in posizione di pivot l'elemento di modulo massimo della colonna) dà $PA=LU$ con $|l_{ik}|\le1$ ed esiste per ogni $A$ invertibile. Se $A$ è simmetrica definita positiva esiste $A=LL^T$ (Cholesky), con $l_{jj}=\sqrt{a_{jj}-\sum_{k<j}l_{jk}^2}$ e $l_{ij}=\frac{a_{ij}-\sum_{k<j}l_{ik}l_{jk}}{l_{jj}}$; costo $\frac13n^3$, senza pivoting.Fattorizzazione LU, pivoting e Cholesky →), cioè meno di un secondo a $10^9$ flop/s.

Errore assoluto e relativo

Sia $x$ il valore esatto e $\tilde x$ un'approssimazione.

Definizione (errori).

• Errore assoluto: $E_a = |x - \tilde x|$ (ha le unità di $x$).
• Errore relativo: $E_r = \dfrac{|x - \tilde x|}{|x|}$ per $x \ne 0$ (adimensionale).

L'errore relativo è quello che conta: sbagliare di $1\,\text{m}$ nella lunghezza di un ponte di $100\,\text{m}$ è poco ($1\%$), sbagliare di $1\,\text{m}$ nella lunghezza di una vite da $5\,\text{cm}$ è assurdo. Un errore relativo di circa $5 \cdot 10^{-k}$ corrisponde a circa $k$ cifre significative corrette.

Esempio. $x = \pi = 3{,}14159265$ e $\tilde x = 3{,}14$. $E_a = 0{,}00159265$, $E_r = 0{,}00159265 / 3{,}14159265 \approx 5{,}07 \cdot 10^{-4}$: circa $3$ cifre corrette.

L'insieme dei numeri di macchina

Definizione (sistema di numeri di macchina $F(\beta, t, L, U)$). Un numero normalizzato è $$x = \pm\,(d_0{,}d_1 d_2 \dots d_{t-1})_\beta \cdot \beta^{\,e}, \qquad d_0 \ne 0,\quad L \le e \le U$$ dove $\beta$ è la base, $t$ il numero di cifre significative (la precisione), $e$ l'esponente. Si aggiunge lo $0$.

In doppia precisione IEEE: $\beta = 2$, $t = 53$ (52 bit di mantissa più l'$1$ nascosto), $e \in [-1022, 1023]$. In singola: $t = 24$, $e \in [-126, 127]$.

Quanti numeri? Segno: $2$ scelte; prima cifra: $\beta - 1$ scelte; le altre $t - 1$ cifre: $\beta^{t-1}$ scelte; esponente: $U - L + 1$ valori; più lo zero: $$\#F = 2(\beta - 1)\beta^{\,t-1}(U - L + 1) + 1.$$

Esempio (sistema "giocattolo"). $F(2, 3, -1, 2)$: numeri $\pm(1{,}m_1 m_2)_2 \cdot 2^e$ con $e \in \{-1, 0, 1, 2\}$. I valori positivi sono $(1 + m/4)\cdot 2^e$ con $m = 0, 1, 2, 3$:

Sono $16$ numeri positivi: $\#F = 2\cdot 1\cdot 4 \cdot 4 + 1 = 33$ con gli opposti e lo zero ✓. Il più grande è $x_{\max} = 7$, il più piccolo positivo $x_{\min} = 0{,}5$.

Nell'intervallo tra $0$ e $x_{\min}$ non c'è nessun numero (nei sistemi reali lo riempiono i denormalizzati: vedi Numeri in virgola mobile IEEE 754Virgola fissa e virgola mobile; formato IEEE 754 a 32 e 64 bit (segno, esponente in eccesso 127, mantissa con 1 implicito); conversioni svolte nei due versi; valori speciali, denormalizzati, intervallo, precisione e arrotondamento.Numeri in virgola mobile IEEE 754 →).

Arrotondamento e precisione di macchina

Un reale $x$ con $x_{\min} \le |x| \le x_{\max}$ non è in generale un numero di macchina: lo si sostituisce con il numero di macchina più vicino, che indichiamo $\mathrm{fl}(x)$ (arrotondamento).

Teorema (errore di arrotondamento). Se $x_{\min} \le |x| \le x_{\max}$, allora $$\mathrm{fl}(x) = x(1 + \delta), \qquad |\delta| \le u = \tfrac12\,\beta^{\,1-t}.$$ $u$ si chiama unità di arrotondamento (unit roundoff); $\varepsilon_M = \beta^{1-t} = 2u$ è la precisione di macchina (distanza tra $1$ e il numero di macchina successivo).

Dimostrazione. Sia $x = (d_0{,}d_1 \dots)_\beta \cdot \beta^e$. In quell'intervallo di esponente i numeri di macchina sono equispaziati con passo $\beta^{\,e+1-t}$ (il peso dell'ultima cifra). Il più vicino dista al massimo mezzo passo: $|x - \mathrm{fl}(x)| \le \tfrac12 \beta^{\,e+1-t}$. Poiché $d_0 \ne 0$ si ha $|x| \ge \beta^{e}$, quindi $$\frac{|x - \mathrm{fl}(x)|}{|x|} \le \frac{\tfrac12\beta^{\,e+1-t}}{\beta^{e}} = \tfrac12\beta^{\,1-t}. \qquad \blacksquare$$

Il punto chiave: l'errore relativo è limitato da $u$, non l'errore assoluto. L'errore assoluto cresce con $|x|$ (nel sistema giocattolo, il passo è $0{,}125$ vicino a $0{,}5$ ma $1$ vicino a $7$).

Esempio (arrotondamento nel sistema giocattolo). Qui $\varepsilon_M = 2^{1-3} = 0{,}25$ e $u = 0{,}125$.

• $x = 2{,}7$ sta tra $2{,}5$ e $3$: distanze $0{,}2$ e $0{,}3$, quindi $\mathrm{fl}(2{,}7) = 2{,}5$. Errore relativo $0{,}2/2{,}7 \approx 0{,}074 \le 0{,}125$ ✓.
• $x = 5{,}4$ sta tra $5$ e $6$: $\mathrm{fl}(5{,}4) = 5$, errore relativo $0{,}4/5{,}4 \approx 0{,}074 \le 0{,}125$ ✓.
• $x = 7{,}6 > 7 + \tfrac12$: overflow (risultato $\infty$ o errore). $x = 0{,}1 < x_{\min}$: underflow (si arrotonda a $0$ o al numero più vicino).
• $1 + 0{,}1 = 1{,}1$ sta tra $1$ e $1{,}25$, più vicino a $1$: $\mathrm{fl}(1{,}1) = 1$. Il $0{,}1$ è assorbito: perché il risultato cambi, a $1$ si deve sommare più di $\varepsilon_M/2=u$ (con $u$ esatto si cade a metà tra $1$ e $1{,}25$ e il pareggio va verso $1$, il numero con ultima cifra pari).

Aritmetica di macchina: le operazioni

Lo standard IEEE richiede che ogni operazione elementare $\circ \in \{+, -, \times, /\}$ tra numeri di macchina dia il risultato esatto arrotondato:

Modello standard. $\mathrm{fl}(a \circ b) = (a \circ b)(1 + \delta)$, con $|\delta| \le u$ (purché non ci sia overflow o underflow).

Ogni operazione introduce quindi un errore relativo piccolo, ma le proprietà algebriche si perdono:

• La somma non è associativa: in doppia precisione $(10^{16} + 1) - 10^{16} = 0$ ma $(10^{16} - 10^{16}) + 1 = 1$, perché vicino a $10^{16}$ il passo tra numeri di macchina è $2$ e il $+1$ viene assorbito.
• $0{,}1 + 0{,}2 \ne 0{,}3$: $0{,}1$ ha sviluppo binario periodico e viene arrotondato (risultato $0{,}30000000000000004$).
• Quindi non si confrontano mai due reali con ==: si controlla $|a - b| \le \text{tol}\,\max(|a|, |b|)$.

Esempio (l'ordine conta, in $F(10,4)$ con arrotondamento a $4$ cifre significative). Si vuole sommare $1234 + 0{,}4 + 0{,}4$. Il risultato esatto è $1234{,}8$.

• Da sinistra: $1234 + 0{,}4 = 1234{,}4 \to 1234$; poi $1234 + 0{,}4 \to 1234$. Risultato $1234$ (errore $0{,}8$).
• Sommando prima i piccoli: $0{,}4 + 0{,}4 = 0{,}8$; poi $1234 + 0{,}8 = 1234{,}8 \to 1235$. Risultato $1235$ (errore $0{,}2$).

Regola pratica: si sommano prima i termini di modulo piccolo.

Cancellazione numerica

Il caso più insidioso è la sottrazione di due numeri quasi uguali. Non è la sottrazione a introdurre errore (se $\frac b2 \le a \le 2b$ la differenza di due numeri di macchina è calcolata esattamente): è che mette in evidenza gli errori già presenti nei dati.

Siano $\tilde a = a(1 + \varepsilon_a)$ e $\tilde b = b(1 + \varepsilon_b)$ le due quantità con i loro errori relativi (ad esempio di arrotondamento). Allora $$\tilde a - \tilde b - (a - b) = a\varepsilon_a - b\varepsilon_b \quad\Longrightarrow\quad \frac{|\tilde a - \tilde b - (a - b)|}{|a - b|} \le \frac{|a| + |b|}{|a - b|}\,\max(|\varepsilon_a|, |\varepsilon_b|).$$ Se $a \approx b$ il fattore $\frac{|a| + |b|}{|a - b|}$ è enorme: l'errore relativo del risultato è molto più grande di quello dei dati. Le cifre in comune ai due numeri si cancellano, e restano solo le ultime, quelle già inquinate.

Procedura (evitare la cancellazione). Se un'espressione contiene la differenza di quantità quasi uguali, la si riscrive in forma algebricamente equivalente senza sottrazione (razionalizzazione, identità trigonometriche, sviluppo di Taylor, formula alternativa).

Esempio 1: differenza di radici. Si calcola $\sqrt{1001} - \sqrt{1000}$ lavorando con $5$ cifre significative.

• $\sqrt{1001} = 31{,}6386\ldots \to 31{,}639$ e $\sqrt{1000} = 31{,}6228\ldots \to 31{,}623$. Ciascuno ha errore relativo piccolo, al più $u = 5 \cdot 10^{-5}$.
• Differenza: $31{,}639 - 31{,}623 = 0{,}016$, con una sola cifra significativa. Il valore esatto è $0{,}0158074\ldots$: errore relativo $\approx 1{,}2 \cdot 10^{-2}$, cioè centinaia di volte l'errore dei dati. Il fattore massimo previsto è $\frac{31{,}639 + 31{,}623}{0{,}016} \approx 3950$.
• Riformulando: $\sqrt{1001} - \sqrt{1000} = \dfrac{1001 - 1000}{\sqrt{1001} + \sqrt{1000}} = \dfrac{1}{31{,}639 + 31{,}623} = \dfrac{1}{63{,}262} = 0{,}015807$: corretto fino all'ultima cifra. La somma al denominatore non ha cancellazione.

Esempio 2: $\dfrac{1 - \cos x}{x^2}$ per $x$ piccolo. Per $x \to 0$ il limite è $\tfrac12$ (Taylor: $\cos x = 1 - \frac{x^2}2 + \dots$, vedi Sviluppi di Mac-Laurin notevoliTabella degli sviluppi di Mac-Laurin da sapere a memoria (e^x, sin, cos, log(1+x), (1+x)^alpha, arctan, sinh, cosh, tan) e regole per combinarli: algebra degli o piccoli, prodotti, funzioni composte, quanti termini tenere.Sviluppi di Mac-Laurin notevoli →). Ma già per $x = 10^{-6}$ il coseno vale $1 - 5 \cdot 10^{-13}$ e nel calcolo si perdono 12 cifre su 16; per $x = 10^{-8}$ il coseno è arrotondato a esattamente $1$ e la formula dà $0$. Si usa invece l'identità $1 - \cos x = 2\sin^2\frac x2$, senza sottrazione: $$\frac{1 - \cos x}{x^2} = \frac{2\sin^2(x/2)}{x^2}.$$

Esempio 3: equazione di secondo grado. Per $ax^2 + bx + c = 0$ con $b^2 \gg 4ac$ e $b > 0$, la radice $x_2 = \frac{-b + \sqrt{b^2 - 4ac}}{2a}$ somma due numeri quasi opposti. Si calcola prima la radice senza cancellazione, $x_1 = \frac{-b - \sqrt{b^2 - 4ac}}{2a}$, e poi si ricava l'altra dal prodotto delle radici $x_1 x_2 = \frac ca$: $$x_2 = \frac{c}{a\,x_1}.$$

Codice

import numpy as np

print(np.finfo(float).eps)                    # precisione di macchina
print(np.spacing(1e8), np.spacing(1e16))      # passo tra numeri di macchina
print(0.1 + 0.2 == 0.3)
print((1e16 + 1) - 1e16, (1e16 - 1e16) + 1)

for k in (2, 4, 6, 8):                        # (1 - cos x)/x^2: formula diretta e riformulata
    x = 10.0**-k
    print(x, (1 - np.cos(x)) / x**2, 2*np.sin(x/2)**2 / x**2)

b, c = 1e8, 1.0                               # x^2 - b x + c = 0, radici ~1e8 e ~1e-8
d = np.sqrt(b*b - 4*c)
x1 = (b + d) / 2
print(x1, (b - d) / 2, c / x1)                # radice grande, radice piccola (cancellazione), radice piccola stabile

Output:

2.220446049250313e-16
1.4901161193847656e-08 2.0
False
0.0 1.0
0.01 0.4999958333473664 0.49999583334722214
0.0001 0.4999999969612645 0.4999999995833334
1e-06 0.5000444502911705 0.4999999999999583
1e-08 0.0 0.5
100000000.0 7.450580596923828e-09 1e-08

La radice piccola con la formula diretta vale $7{,}45 \cdot 10^{-9}$ invece di $10^{-8}$ (errore del $25\%$); con $x_2 = c/(a x_1)$ è esatta. In MATLAB: eps (precisione di macchina), eps(x) (passo vicino a x), realmax, realmin, sqrt, cos.

Errori tipici

• Confondere $u$ ed $\varepsilon_M$ (differiscono per un fattore $2$): $\varepsilon_M$ è la distanza tra $1$ e il successivo, $u$ è l'errore relativo massimo di arrotondamento.
• Dire che l'errore assoluto di arrotondamento è sempre $\le u$: è l'errore relativo a essere limitato.
• Dare la colpa alla sottrazione: la sottrazione è esatta, è la perdita di cifre significative nei dati a diventare visibile.
• Confrontare reali con == o fermare un ciclo quando due numeri sono uguali.

Domande d'esame

1. Definire errore assoluto e relativo e dire quale dei due è limitato dall'arrotondamento. Risposta: il relativo, $|\delta| \le u = \frac12\beta^{1-t}$.
2. Quanto vale la precisione di macchina in doppia precisione e cosa rappresenta? Risposta: $\varepsilon_M = 2^{-52} \approx 2{,}22\cdot 10^{-16}$: distanza tra $1$ e il numero successivo; l'errore relativo di arrotondamento è al più $u = \varepsilon_M/2$.
3. Perché $\sqrt{x+1} - \sqrt{x}$ è calcolata male per $x$ grande e come si migliora? Risposta: cancellazione tra numeri quasi uguali; si usa $\frac{1}{\sqrt{x+1} + \sqrt{x}}$.
4. Quali interi si rappresentano con $16$ bit in complemento a due e cosa succede sommando $1$ al massimo? Risposta: $[-32768,\,32767]$; $32767+1$ dà overflow e il risultato si avvolge a $-32768$ (in MATLAB int16 satura a $32767$).
5. L'addizione di macchina è associativa? Risposta: no, per l'assorbimento: $(10^{16} + 1) - 10^{16} = 0$ ma $(10^{16} - 10^{16}) + 1 = 1$.