12 - NumPy: array, vettorizzazione e broadcasting

Corso Machine Learning Engineer · Modulo 01 · Python per i dati · Lezione 5

Hai un ciclo che somma i quadrati di un milione di numeri e impiega un decimo di secondo, e un'espressione di una riga che fa lo stesso lavoro in meno di un millesimo. Questa lezione, che viene dopo i numeri in virgola mobile, spiega da dove nasce quella differenza e ti dà le tre idee che reggono quasi tutto il codice di machine learning: le viste, che non copiano, il broadcasting, che allinea forme diverse, e gli assi, che dicono su quale dimensione si calcola. Alla fine sai prevedere la forma del risultato di un'operazione prima di eseguirla e riconosci l'errore di forma che produce un numero sbagliato senza alcun messaggio.

Il luogo comune dice che NumPy è veloce perché è scritto in C. È solo metà della verità: una funzione C chiamata un milione di volte da un ciclo Python resta lenta. NumPy è veloce perché i dati sono omogenei e contigui in memoria e perché un'unica chiamata lavora su tutti gli elementi. Conviene impararlo bene perché quasi ogni libreria di dati (Pandas, scikit-learn, Matplotlib) scambia numeri sotto forma di array NumPy, o di oggetti che gli somigliano.

Perché una lista Python non basta

L'analogia: la cassetta delle lettere contro la scatola delle uova

Una lista Python è una fila di biglietti. Ogni biglietto non contiene il numero: contiene l'indirizzo di una scatola in magazzino, dove sta un oggetto float completo (con il suo contatore di riferimenti, il suo tipo, il suo valore). Per sommare mille numeri l'interprete deve andare mille volte a cercare una scatola in giro per il magazzino, aprirla, controllare di che tipo è.

Un array NumPy è una scatola di uova: tutte dello stesso tipo, una accanto all'altra, in un blocco unico. Per sommarle basta scorrere il blocco. Non c'è nulla da controllare, perché la scatola ha un'etichetta sola sul coperchio: «uova da 8 byte».

Misuriamolo

import sys
import numpy as np

lista = [float(i) for i in range(1_000_000)]
array = np.arange(1_000_000, dtype=np.float64)

byte_lista = sys.getsizeof(lista) + sum(sys.getsizeof(x) for x in lista)
print("lista: byte dei puntatori      :", sys.getsizeof(lista))
print("lista: byte di un singolo float:", sys.getsizeof(lista[1]))
print("lista: totale con gli oggetti  :", byte_lista)
print("array: byte                    :", array.nbytes)
print("rapporto di memoria            :", round(byte_lista / array.nbytes, 1))

Output ottenuto:

lista: byte dei puntatori      : 8448728
lista: byte di un singolo float: 24
lista: totale con gli oggetti  : 32448728
array: byte                    : 8000000
rapporto di memoria            : 4.1

La lista spende 8 byte per ogni puntatore e altri 24 per ogni oggetto float, quindi circa 32 byte per numero contro 8. Oltre allo spazio, conta la località: un blocco contiguo si legge in sequenza, e la memoria del computer è molto più veloce a leggere in sequenza che a saltare da un indirizzo all'altro.

Trade-off. Un array ha un solo tipo e una dimensione fissa. Se devi mescolare stringhe, numeri e None, o aggiungere elementi di continuo, una lista (o un DataFrame) è lo strumento giusto. NumPy paga la sua velocità con la rigidità.

Anatomia di un array: shape, dtype, strides

Un ndarray è, concretamente, un blocco di memoria più tre numeri che dicono come leggerlo:

• dtype: che tipo di numero è ogni cella (float64, int32, bool...). Fissa anche quanti byte occupa (itemsize).
• shape: quante celle ci sono lungo ogni dimensione, per esempio (3, 4) per 3 righe e 4 colonne.
• strides: di quanti byte ci si sposta in memoria per passare all'elemento successivo lungo ciascuna dimensione.

I primi due numeri dicono che cosa c'è nel blocco, il terzo dice come attraversarlo. Separare le due cose è la scelta di progetto che rende economiche trasposte, fette e cambi di forma: nessuno di questi gesti tocca i dati, cambia solo il modo di leggerli. Il legame fra i tre è una formula sola.

La formula che lega tutto è l'indirizzo dell'elemento a[i, j]:

indirizzo(i, j) = base + i · strides[0] + j · strides[1]

Per una matrice (3, 4) di float64 memorizzata per righe, passare alla colonna successiva sposta di 8 byte (un float64), passare alla riga successiva sposta di 4 · 8 = 32 byte. Quindi strides = (32, 8).

import numpy as np

x = np.arange(12, dtype=np.float64).reshape(3, 4)
print(x)
print("shape:", x.shape, "| dtype:", x.dtype, "| itemsize:", x.itemsize, "| strides:", x.strides)

t = x.T                      # la trasposta
print("trasposta strides:", t.strides, "| condivide la memoria con x:", np.shares_memory(x, t))

Output ottenuto:

[[ 0.  1.  2.  3.]
 [ 4.  5.  6.  7.]
 [ 8.  9. 10. 11.]]
shape: (3, 4) | dtype: float64 | itemsize: 8 | strides: (32, 8)
trasposta strides: (8, 32) | condivide la memoria con x: True

La trasposta non sposta nemmeno un numero: scambia i due strides da (32, 8) a (8, 32). Cambia il modo di leggere, non il contenuto. Questo è il motivo per cui x.T costa praticamente zero, e anche il motivo per cui NumPy può permettersi tante operazioni «gratuite» di cambio di forma.

reshape funziona allo stesso modo: se il numero totale di elementi torna, ricalcola shape e strides senza toccare i dati (quando la disposizione in memoria lo permette).

Creare array

import numpy as np

print(np.zeros((2, 3)))
print(np.ones(3, dtype=np.int64))
print(np.eye(3))
print(np.arange(0, 1, 0.25))          # passo fisso, estremo destro escluso
print(np.linspace(0, 1, 5))           # numero di punti fisso, estremi inclusi
print(np.full((2, 2), 7.5))

Output ottenuto:

[[0. 0. 0.]
 [0. 0. 0.]]
[1 1 1]
[[1. 0. 0.]
 [0. 1. 0.]
 [0. 0. 1.]]
[0.   0.25 0.5  0.75]
[0.   0.25 0.5  0.75 1.  ]
[[7.5 7.5]
 [7.5 7.5]]

Trappola reale: arange con passo decimale. np.arange(inizio, fine, passo) produce ceil((fine − inizio) / passo) elementi (è la formula della documentazione di numpy.arange, che avverte anche che l'ultimo elemento può risultare maggiore di stop). Se il quoziente è un intero «per un soffio» in virgola mobile, l'esito sorprende:

import numpy as np

a = np.arange(0, 2.1, 0.3)
print("elementi:", len(a), "| ultimo:", a[-1], "| quoziente 2.1/0.3 =", 2.1 / 0.3)
print("linspace, 8 punti esatti:", np.linspace(0, 2.1, 8))

Output ottenuto:

elementi: 8 | ultimo: 2.1 | quoziente 2.1/0.3 = 7.000000000000001
linspace, 8 punti esatti: [0.  0.3 0.6 0.9 1.2 1.5 1.8 2.1]

Sulla carta 2.1 / 0.3 = 7 esatto, quindi 7 elementi (da 0 a 1.8) e l'estremo destro escluso. Ma 0.3 e 2.1 non sono rappresentabili (lezione precedente): il quoziente calcolato vale 7.000000000000001, il soffitto è 8, e arange restituisce otto elementi, con l'ultimo che coincide con l'estremo che doveva essere escluso. Quando conosci quanti punti vuoi, usa np.linspace(inizio, fine, n): il numero di punti è un parametro, non un risultato dell'aritmetica in virgola mobile.

Tipi numerici e overflow silenzioso

Ogni dtype intero ha un intervallo finito. Con int8 (8 bit con segno) i valori vanno da −128 a 127. Cosa succede oltre? Il contatore si riavvolge, come il contachilometri di un'auto vecchia.

import numpy as np

piccolo = np.array([127], dtype=np.int8)
print("127 + 1 in int8 :", piccolo + 1)
print("intervallo int8 :", np.iinfo(np.int8).min, np.iinfo(np.int8).max)

img_a = np.array([200], dtype=np.uint8)      # come il valore di un pixel
img_b = np.array([100], dtype=np.uint8)
print("uint8: 200 + 100          :", img_a + img_b, "(riavvolto: 300 - 256)")
print("uint8 convertito prima    :", img_a.astype(np.int64) + img_b.astype(np.int64))

print("float32 di 0.1  :", float(np.float32(0.1)))
print("float64 di 0.1  :", float(np.float64(0.1)))

Output ottenuto:

127 + 1 in int8 : [-128]
intervallo int8 : -128 127
uint8: 200 + 100          : [44] (riavvolto: 300 - 256)
uint8 convertito prima    : [300]
float32 di 0.1  : 0.10000000149011612
float64 di 0.1  : 0.1

Sugli array, nessun errore e nessun avviso: solo un numero sbagliato. (Sui singoli scalari NumPy può emettere un RuntimeWarning, ma i calcoli su array, quelli che contano nei dati, si riavvolgono in silenzio.) La derivazione è semplice: con n bit si hanno 2ⁿ valori; in aritmetica modulo 2ⁿ, 300 mod 256 = 44. Le immagini a 8 bit sono il caso classico: sommare due immagini uint8 senza convertirle prima produce pixel scuri dove dovevano essere chiari.

Il float32 occupa metà della memoria del float64 ma ha circa 7 cifre decimali significative contro 15-16. Per le reti neurali spesso basta (e dimezza memoria e tempo); per sommare milioni di termini o per calcolare varianze spesso no, come visto nella lezione sui numeri in virgola mobile.

Regola pratica: lascia fare a NumPy il tipo di default (float64, int64) finché non misuri un problema di memoria o di velocità. Ottimizzare il tipo prima di avere un problema aggiunge rischio senza un beneficio dimostrato.

Indicizzazione: viste contro copie

Questa è la sezione che evita più bug. Esistono tre modi di selezionare, con comportamenti diversi:

La regola completa è nella guida Copies and views: l'indicizzazione di base (le fette) restituisce sempre una vista, l'indicizzazione avanzata (liste di indici, maschere) sempre una copia.

Analogia. Una vista è una finestra su una stanza: se qualcuno sposta un mobile, lo vedi cambiare dalla finestra, e se lo sposti tu dalla finestra, cambia nella stanza. Una copia è una fotografia: da quel momento le due vite sono indipendenti.

import numpy as np

a = np.arange(10)
vista = a[2:5]
vista[0] = 99                       # modifico la "finestra"
print("originale dopo la modifica della vista:", a)

copia = a[[2, 3, 4]]
copia[0] = -1                       # modifico la "fotografia"
print("originale dopo la modifica della copia :", a)

print("vista condivide memoria:", np.shares_memory(a, vista))
print("copia condivide memoria:", np.shares_memory(a, copia))

indipendente = a[2:5].copy()        # per avere una copia esplicita
print("copia esplicita condivide:", np.shares_memory(a, indipendente))

Output ottenuto:

originale dopo la modifica della vista: [ 0  1 99  3  4  5  6  7  8  9]
originale dopo la modifica della copia : [ 0  1 99  3  4  5  6  7  8  9]
vista condivide memoria: True
copia condivide memoria: False
copia esplicita condivide: False

Quando serve cosa. La vista è velocissima e non occupa memoria, perfetta per leggere. Se vuoi modificare un sottoinsieme senza toccare l'originale, chiama .copy() esplicitamente. In caso di dubbio, np.shares_memory(a, b) risponde con certezza.

Vettorizzazione: un'istruzione, un milione di elementi

Vettorizzare significa sostituire un ciclo Python con un'unica espressione che NumPy esegue internamente in codice compilato, su tutti gli elementi.

Confrontiamo due modi di calcolare la somma dei quadrati di un milione di numeri.

import numpy as np

import time

v = np.random.default_rng(0).random(1_000_000)

t0 = time.perf_counter()
s_ciclo = 0.0
for e in v.tolist():
    s_ciclo += e * e
t_ciclo = time.perf_counter() - t0

t0 = time.perf_counter()
s_vett = np.dot(v, v)
t_vett = time.perf_counter() - t0

print(f"ciclo Python : {t_ciclo:.4f} s")
print(f"vettorizzato : {t_vett:.5f} s")
print(f"rapporto     : circa {t_ciclo / t_vett:.0f} volte")
print("stesso risultato (entro la tolleranza):", np.isclose(s_ciclo, s_vett))

Output ottenuto:

ciclo Python : 0.1119 s
vettorizzato : 0.00088 s
rapporto     : circa 128 volte
stesso risultato (entro la tolleranza): True

I tempi cambiano da una esecuzione all'altra e da una macchina all'altra, quindi l'uscita sopra è una sola esecuzione fra cinque. In cinque esecuzioni successive su un Ryzen 5 5600H, Python 3.14.4 e NumPy 2.3.5, il ciclo ha impiegato da 0,10 a 0,15 s, la versione vettorizzata circa 0,0009 s, e il rapporto è stato fra 105 e 168 volte: l'ordine di grandezza (un fattore intorno a cento) è stabile, la cifra esatta no. Nota che il ciclo parte già da v.tolist(), cioè con la conversione in oggetti Python fuori dalla misura: è un confronto favorevole al ciclo.

Il motivo è quello dell'analogia iniziale. Per ogni elemento, il ciclo Python deve leggere l'oggetto float, controllarne il tipo, trovare la funzione di moltiplicazione, creare un nuovo oggetto per il risultato e ripetere per la somma. NumPy, con dtype fissato in anticipo, esegue un ciclo compilato che fa solo le moltiplicazioni e le somme.

Le funzioni che operano elemento per elemento e che NumPy chiama ufunc (universal functions) sono ovunque: np.exp, np.log, np.sqrt, +, *, i confronti. Si compongono liberamente:

import numpy as np

x = np.linspace(-3, 3, 7)
sigmoide = 1 / (1 + np.exp(-x))          # definizione: 1 / (1 + e^{-x})
print(np.round(x, 2))
print(np.round(sigmoide, 4))

Output ottenuto:

[-3. -2. -1.  0.  1.  2.  3.]
[0.0474 0.1192 0.2689 0.5    0.7311 0.8808 0.9526]

Una sola riga, nessun ciclo, funziona per 7 elementi come per 7 milioni.

Quando NON vettorizzare. Se il calcolo di ogni elemento dipende dal precedente (una ricorrenza), non è parallelizzabile in blocco e un ciclo può essere inevitabile (esistono strumenti come Numba per compilarlo, ma sono fuori da questa lezione). E una vettorizzazione che crea array intermedi enormi può esaurire la memoria: un ciclo a blocchi è a volte la risposta giusta.

Assi e aggregazioni

Una funzione di aggregazione (sum, mean, max, std...) riduce un array. Su quale dimensione? Lo decide axis. La regola da ricordare è una sola:

L'asse che indichi è quello che scompare.

Per una matrice con 3 righe e 4 colonne, axis=0 collassa le righe (ottieni un valore per colonna: forma (4,)), axis=1 collassa le colonne (un valore per riga: forma (3,)).

import numpy as np

m = np.array([[1, 2, 3, 4],
              [5, 6, 7, 8],
              [9, 10, 11, 12]])
print("somma di tutto     :", m.sum())
print("axis=0 (per colonna):", m.sum(axis=0), "forma", m.sum(axis=0).shape)
print("axis=1 (per riga)   :", m.sum(axis=1), "forma", m.sum(axis=1).shape)
print("keepdims            :", m.sum(axis=1, keepdims=True).shape)

Output ottenuto:

somma di tutto     : 78
axis=0 (per colonna): [15 18 21 24] forma (4,)
axis=1 (per riga)   : [10 26 42] forma (3,)
keepdims            : (3, 1)

keepdims=True conserva l'asse con lunghezza 1: serve spesso per il broadcasting che vediamo tra poco. Nei dati tabulari, la convenzione è righe = esempi, colonne = caratteristiche: X.mean(axis=0) è la media di ciascuna caratteristica.

Broadcasting: le regole, derivate

L'analogia: la stampante che ripete l'etichetta

Hai 100 cartelle e un'unica etichetta con la data di oggi. Non stampi 100 etichette: la applichi a tutte le cartelle. Il broadcasting fa lo stesso: quando due array hanno forme diverse, NumPy «stira» virtualmente quello più piccolo, ripetendolo, senza copiare i dati.

Le regole

Per combinare due array, NumPy confronta le forme dall'ultima dimensione verso la prima. Due dimensioni sono compatibili se:

1. sono uguali, oppure
2. una delle due vale 1 (che viene allungata fino a combaciare con l'altra).

Se manca una dimensione (un array ha meno dimensioni), si considera di lunghezza 1 a sinistra. Se qualche coppia non rispetta nessuna delle due regole: errore (ValueError, come descrive la guida Broadcasting). La forma del risultato è, in ogni posizione, il valore più grande fra i due, con una sola eccezione: una dimensione di lunghezza 0 contro una di lunghezza 1 dà 0 (per esempio (3, 1) e (0,) danno (3, 0), verificato).

Applichiamole a tre casi, leggendo le forme da destra:

import numpy as np

X = np.array([[1., 2., 3., 4.],
              [5., 6., 7., 8.],
              [9., 10., 11., 12.]])
media = X.mean(axis=0)                  # forma (4,)
print("X - media:\n", X - media)        # (3,4) - (4,) -> (3,4)

colonna = np.array([[10], [20], [30]])  # (3,1)
riga = np.array([[1, 2, 3, 4]])         # (1,4)
print("tabella di somme:\n", colonna + riga)

v = np.array([1., 2., 3.])              # (3,)
try:
    X - v
except ValueError as e:
    print("errore atteso:", e)
print("sistemato con v[:, None]:\n", X - v[:, None])   # (3,4) - (3,1)

Output ottenuto:

X - media:
 [[-4. -4. -4. -4.]
 [ 0.  0.  0.  0.]
 [ 4.  4.  4.  4.]]
tabella di somme:
 [[11 12 13 14]
 [21 22 23 24]
 [31 32 33 34]]
errore atteso: operands could not be broadcast together with shapes (3,4) (3,) 
sistemato con v[:, None]:
 [[0. 1. 2. 3.]
 [3. 4. 5. 6.]
 [6. 7. 8. 9.]]

Il trucco v[:, None] (equivalente a v.reshape(3, 1)) inserisce un asse di lunghezza 1 e trasforma un vettore riga in colonna. È il gesto più frequente del broadcasting.

Perché non costa memoria

Il segreto è nei strides: lungo un asse «stirato», lo stride vale zero. Il passo in memoria è nullo, quindi si rilegge sempre lo stesso elemento.

import numpy as np

stirato = np.broadcast_to(np.arange(4), (3, 4))
print("forma:", stirato.shape, "| strides:", stirato.strides, "(zero byte per passare alla riga successiva)")

Output ottenuto:

forma: (3, 4) | strides: (0, 8) (zero byte per passare alla riga successiva)

Standardizzare una matrice

Un'operazione che incontrerai in quasi tutti i modelli è la standardizzazione: sottrarre a ogni colonna la sua media e dividere per la sua deviazione standard, z = (x − μ) / σ. Con il broadcasting è una riga:

import numpy as np

rng = np.random.default_rng(1)
dati = rng.normal(loc=[50, 1000, 0.5], scale=[10, 300, 0.1], size=(500, 3))

z = (dati - dati.mean(axis=0)) / dati.std(axis=0)
print("media di ogni colonna ~ 0:", np.allclose(z.mean(axis=0), 0))
print("std di ogni colonna   = 1:", np.allclose(z.std(axis=0), 1))

Output ottenuto:

media di ogni colonna ~ 0: True
std di ogni colonna   = 1: True

Le tre colonne partivano da scale diversissime (decine, migliaia, decimi). Dopo la standardizzazione hanno tutte media 0 e deviazione standard 1.

Distanze a coppie senza cicli, e il costo nascosto

Date n punti in d dimensioni, la distanza euclidea tra ogni coppia si ottiene estendendo gli assi:

import numpy as np

rng = np.random.default_rng(1)
punti = rng.normal(size=(4, 2))
diff = punti[:, None, :] - punti[None, :, :]       # (4,1,2) - (1,4,2) -> (4,4,2)
distanze = np.sqrt((diff ** 2).sum(axis=-1))        # somma sull'ultimo asse -> (4,4)
print(np.round(distanze, 3))
print("diagonale (distanza di un punto da sé):", distanze.diagonal())

Output ottenuto:

[[0.    2.125 0.674 0.915]
 [2.125 0.    1.842 2.074]
 [0.674 1.842 0.    1.449]
 [0.915 2.074 1.449 0.   ]]
diagonale (distanza di un punto da sé): [0. 0. 0. 0.]

Trade-off. L'array intermedio diff ha forma (n, n, d): per n = 100 000 punti e d = 10 servirebbero 10¹⁰ · 10 · 8 byte, cioè 800 GB. Il broadcasting è elegante, ma la memoria non è infinita. Per insiemi grandi si usano algoritmi dedicati (per esempio scipy.spatial.distance.cdist, o il calcolo a blocchi).

Prodotto matriciale contro prodotto elemento per elemento

L'operatore * moltiplica elemento per elemento (con broadcasting). Il prodotto matriciale si scrive @ (o np.matmul, np.dot). Per A di forma (m, n) e B di forma (n, p):

(A @ B)ᵢⱼ = Σₖ Aᵢₖ · Bₖⱼ, con risultato di forma (m, p).

Il numero di colonne di A deve coincidere con il numero di righe di B: è la condizione per poter fare la somma su k.

import numpy as np

A = np.array([[1, 2], [3, 4]])
B = np.array([[0, 1], [1, 0]])
print("elemento per elemento A * B:\n", A * B)
print("prodotto matriciale   A @ B:\n", A @ B)
print("B @ A diverso da A @ B:", not np.array_equal(A @ B, B @ A))

Output ottenuto:

elemento per elemento A * B:
 [[0 2]
 [3 0]]
prodotto matriciale   A @ B:
 [[2 1]
 [4 3]]
B @ A diverso da A @ B: True

Confondere * con @ è un errore comunissimo: nessuna eccezione se le forme sono compatibili con entrambi, risultati completamente diversi. Il prodotto matriciale non è commutativo: A @ B e B @ A sono in generale matrici diverse. È l'operazione centrale dell'algebra lineare; qui la usiamo senza svilupparne la teoria.

Maschere, where e valori mancanti

Una maschera booleana è un array di True/False della stessa forma, e seleziona gli elementi veri.

import numpy as np

temperature = np.array([18.5, 22.0, np.nan, 30.1, 15.2, np.nan, 27.8])
print("nan == nan:", np.nan == np.nan)               # NaN non è uguale a nulla, nemmeno a sé stesso
validi = ~np.isnan(temperature)
print("media ingenua      :", temperature.mean())    # il NaN contagia tutto
print("media dei validi   :", temperature[validi].mean())
print("equivalente        :", np.nanmean(temperature))

caldo = np.where(temperature > 25, "caldo", "mite")  # (NaN > 25 è False)
print(caldo)
print("clip tra 20 e 28   :", np.clip(temperature, 20, 28))

Output ottenuto:

nan == nan: False
media ingenua      : nan
media dei validi   : 22.72
equivalente        : 22.72
['mite' 'mite' 'mite' 'caldo' 'mite' 'mite' 'caldo']
clip tra 20 e 28   : [20.  22.   nan 28.  20.   nan 27.8]

NaN («not a number») è il valore che IEEE 754 riserva ai risultati indefiniti. La sua proprietà più insidiosa: qualunque confronto con NaN è falso (compreso NaN == NaN), per cui si testa solo con np.isnan. E si propaga: la media di un array che contiene un solo NaN è NaN. Decidere cosa fare dei valori mancanti (scartarli, riempirli, segnalarli) è una scelta di analisi, non un dettaglio tecnico: ne riparliamo con Pandas.

Numeri casuali riproducibili

Per simulazioni ed esperimenti serve una sorgente di numeri pseudo-casuali riproducibile: con lo stesso seed (seme) si ottiene la stessa sequenza, quindi lo stesso risultato.

import numpy as np

rng = np.random.default_rng(seed=42)       # il modo moderno
print(rng.integers(0, 10, size=5))
print(np.round(rng.normal(size=3), 4))
print(rng.choice(["a", "b", "c"], size=4, replace=True))

rng2 = np.random.default_rng(seed=42)
print("stessa sequenza con lo stesso seed:", rng2.integers(0, 10, size=5))

Output ottenuto:

[0 7 6 4 4]
[ 0.9406 -1.951  -1.3022]
['c' 'c' 'c' 'c']
stessa sequenza con lo stesso seed: [0 7 6 4 4]

Si usa np.random.default_rng (un oggetto Generator) e lo si passa alle funzioni che lo richiedono, invece dello stato globale np.random.seed, che è l'interfaccia più vecchia. Il generatore locale evita che due parti del programma si disturbino a vicenda; la documentazione di NumPy lo raccomanda per il nuovo codice. Con versioni diverse di NumPy la sequenza per lo stesso seed non è garantita identica: la politica di compatibilità del generatore riserva ai progettisti la possibilità di cambiare algoritmo, e perfino chiamare rng.random() cinque volte può non dare gli stessi numeri di rng.random(5). Per questo in questo corso ogni esperimento dichiara le versioni usate.

Errori tipici (sintomo, causa, rimedio)

Il broadcasting silenzioso

È la trappola più costosa, perché non produce nessun errore. Un vettore di previsioni di forma (n, 1) e uno di valori veri di forma (n,) non danno errore: danno una matrice (n, n), e la media di quella matrice è un numero qualunque.

import numpy as np

y_vero = np.array([1.0, 2.0, 3.0])
y_prev = np.array([[1.0], [2.0], [4.0]])          # forma (3,1) invece di (3,)

mse_sbagliato = ((y_prev - y_vero) ** 2).mean()
mse_giusto = ((y_prev[:, 0] - y_vero) ** 2).mean()
print("(3,1) - (3,) ha forma:", (y_prev - y_vero).shape)
print("MSE sbagliato:", round(float(mse_sbagliato), 4), "| MSE giusto:", round(float(mse_giusto), 4))

Output ottenuto:

(3,1) - (3,) ha forma: (3, 3)
MSE sbagliato: 2.3333 | MSE giusto: 0.3333

Rimedio: controlla le forme con assert y_prev.shape == y_vero.shape all'inizio delle funzioni di errore. Una riga che costa nulla.

Le altre trappole

• Modifico una vista senza volerlo. Sintomo: i dati originali cambiano. Rimedio: .copy() quando serve indipendenza.
• Asse sbagliato in un'aggregazione. Sintomo: forma del risultato inattesa. Rimedio: stampa la forma, ricorda «l'asse indicato scompare».
• Overflow di interi piccoli (uint8, int16). Rimedio: converti prima di fare i conti (astype).
• == tra float. Rimedio: np.isclose, np.allclose (vedi la lezione sui numeri in virgola mobile).
• Un ciclo Python su un array. Sintomo: lentezza. Rimedio: cerca l'espressione vettorizzata; se non esiste, valuta se il ciclo è davvero inevitabile.

Le prime due trappole producono un risultato di forma o contenuto inatteso e si scoprono stampando shape; l'overflow e il confronto fra float producono invece numeri plausibili e sbagliati, e vanno prevenuti per costruzione. Il ciclo su un array è l'unica che non dà risposte sbagliate, solo lente: per questo è l'ultima da correggere.

Best practice

• Stampa shape e dtype quando qualcosa non torna: è il primo debug, e risolve la maggioranza dei casi.
• Verifica le forme con assert nelle funzioni che le danno per scontate.
• Preferisci operazioni su array ai cicli, ma non sacrificare la leggibilità: una riga incomprensibile non è un guadagno.
• Usa keepdims=True o [:, None] esplicitamente per rendere visibile il broadcasting.
• Fissa il seed e passa il Generator come argomento.
• Mantieni float64 finché non misuri un problema.

Le sei regole servono lo stesso scopo: rendere visibile ciò che NumPy fa in silenzio. Forme e tipi stampati o verificati con assert trasformano un numero sbagliato in un errore con un nome, e il seed fissato trasforma un esperimento irripetibile in uno che un collega può rifare. Costano poche righe all'inizio di una funzione e risparmiano ore di ricerca.

Mini-progetto

Scrivi una funzione standardizza(X) che:

1. controlla che X abbia due dimensioni (altrimenti solleva ValueError con un messaggio chiaro);
2. restituisce (Z, media, dev_std), dove Z è la matrice standardizzata;
3. gestisce le colonne con deviazione standard zero (colonne costanti) senza produrre NaN: scegli e documenta cosa restituire.

I tre requisiti coprono il controllo dell'ingresso, il contratto di uscita e il caso limite, che è dove le funzioni numeriche di solito si rompono. Verifica con un X di forma (100, 5) in cui una colonna è costante, e controlla che Z.mean(axis=0) sia vicino a zero con np.allclose. Prova prima a scrivere la funzione, poi confrontala con la soluzione commentata qui sotto.

La soluzione commentata del mini-progetto: standardizza

Mettiamo insieme viste, assi e broadcasting nella funzione chiesta dal mini-progetto, con attenzione ai casi limite.

import numpy as np

def standardizza(X):
    """Restituisce (Z, media, dev_std). Le colonne costanti restano a zero invece di diventare NaN."""
    X = np.asarray(X, dtype=float)
    if X.ndim != 2:
        raise ValueError(f"serve una matrice con 2 dimensioni, ricevute {X.ndim}")
    media = X.mean(axis=0)
    dev = X.std(axis=0)
    dev_sicura = np.where(dev == 0, 1.0, dev)           # evita la divisione per zero nelle colonne costanti
    return (X - media) / dev_sicura, media, dev

rng = np.random.default_rng(3)
X = rng.normal(loc=[10.0, -5.0, 3.0, 100.0, 0.0], scale=[2.0, 0.5, 1.0, 30.0, 1.0], size=(100, 5))
X[:, 2] = 7.0                                            # una colonna costante

Z, media, dev = standardizza(X)
print("forma di Z:", Z.shape, "| nessun NaN:", not np.isnan(Z).any())
print("media delle colonne ~ 0:", np.allclose(Z.mean(axis=0), 0))
print("dev. std (colonne non costanti) = 1:", np.allclose(Z.std(axis=0)[[0, 1, 3, 4]], 1))
print("colonna costante -> tutta a 0:", np.allclose(Z[:, 2], 0), "| dev. std originale:", dev[2])

try:
    standardizza(np.array([1.0, 2.0, 3.0]))
except ValueError as e:
    print("ValueError:", e)

# la trasformazione si inverte: Z * dev + media ritrova X (dove la dev. std non era zero)
X_ricostruita = Z * np.where(dev == 0, 1.0, dev) + media
print("X ricostruita da Z:", np.allclose(X_ricostruita, X))

Output ottenuto:

forma di Z: (100, 5) | nessun NaN: True
media delle colonne ~ 0: True
dev. std (colonne non costanti) = 1: True
colonna costante -> tutta a 0: True | dev. std originale: 0.0
ValueError: serve una matrice con 2 dimensioni, ricevute 1
X ricostruita da Z: True

Come leggere la soluzione. np.asarray(X, dtype=float) accetta anche liste e converte in array di decimali senza copiare se è già un array adatto. Il controllo X.ndim != 2 rifiuta subito un vettore: meglio un errore chiaro che un risultato sbagliato per broadcasting. X.mean(axis=0) e X.std(axis=0) hanno forma (5,), e (X - media) / dev_sicura sfrutta il broadcasting di vettori (5,) su una matrice (100, 5).

Il punto delicato è la colonna costante: la sua deviazione standard è zero e 0/0 darebbe NaN in tutta la colonna (con un avviso). np.where(dev == 0, 1.0, dev) sostituisce gli zeri con 1 solo per la divisione: la colonna costante diventa tutta zeri, che è la scelta onesta («non porta informazione»), e dev originale viene restituita intatta così chi chiama sa che quella colonna era costante.

L'ultima verifica mostra che la trasformazione è reversibile: salvare media e dev permette di riportare nuovi dati alla stessa scala o di tornare alle unità originali (è proprio ciò che fa uno scalatore di scikit-learn, dove media e dev si calcolano sul solo training e poi si applicano, immutati, ai dati di test). Una nota di stile: la funzione non modifica X (costruisce un array nuovo) e restituisce tutto ciò che serve per invertire l'operazione; sono le due abitudini che rendono prevedibile il codice sui dati.

Challenge

• Easy. Calcola la media mobile a finestra 3 di un array di 10 numeri, senza cicli, usando solo slicing e somme.
• Medium. Date due matrici A (n, d) e B (m, d), calcola con il broadcasting la matrice (n, m) della distanza al quadrato. Poi dimostra l'identità ‖a − b‖² = ‖a‖² − 2 a·b + ‖b‖² e riscrivila usando un prodotto matriciale per evitare l'array intermedio (n, m, d).
• Hard. Misura con time.perf_counter il tempo della somma dei quadrati per n = 10³, 10⁴, 10⁵, 10⁶, 10⁷ nella versione con ciclo e in quella vettorizzata. Da quale n conviene vettorizzare? Perché per n piccolo il vantaggio si riduce?
• Expert. Con la formula ceil((fine − inizio) / passo), cerca per tentativi una coppia (fine, passo) decimale per cui len(np.arange(0, fine, passo)) differisce dal conteggio esatto (usa fractions.Fraction per il valore esatto) e spiega il perché guardando la rappresentazione binaria del quoziente.
• Enterprise. Prendi una tabella reale della tua organizzazione con almeno cinque colonne numeriche. Standardizzala con NumPy e documenta, per ciascuna colonna, che cosa significa un valore z = 2.

I livelli salgono dalla pratica con le fette (Easy) all'algebra (Medium), alla misura (Hard) e alla ricerca di un controesempio numerico (Expert); l'ultimo non ha una soluzione di codice ma una spiegazione per chi leggerà i tuoi dati. Per Hard ricorda di ripetere ogni misura più volte e di riportare un intervallo, come sopra.

Che cosa sai fare adesso

Sai leggere un array come un blocco di memoria con dtype, shape e strides, e sai dire quando una selezione è una vista e quando una copia. Sai scegliere l'asse di un'aggregazione ricordando che quello indicato scompare, prevedere la forma di un'operazione con il broadcasting allineando le forme da destra, e riconoscere il caso (n, 1) contro (n,) che produce un numero sbagliato senza avvisare. Il passo successivo è Pandas, che mette etichette e tipi misti sopra gli array per lavorare con tabelle vere.

Nella stessa serie trovi anche: 15 - scikit-learn: stimatori, pipeline e validazione, 23 - Algebra lineare 1: vettori, matrici, sistemi, 56 - Reti neurali ricorrenti (RNN): dare una memoria alla rete (che richiamano questa lezione).

Fonti

Riferimento bibliografico non consultato in questa revisione: C. R. Harris et al., Array programming with NumPy, Nature 585, 2020. Tutte le uscite sono state ricalcolate con Python 3.14.4 e NumPy 2.3.5 (ogni blocco eseguito due volte, uscita identica, salvo il blocco dei tempi, riportato come intervallo).