NumPy

Introduzione

NumPy (Numerical Python) è il pacchetto fondamentale per il data science e il calcolo scientifico in Python. Si tratta di una libreria open source che fornisce strutture dati ottimizzate e funzioni per la manipolazione di array multidimensionali, offrendo prestazioni paragonabili a quelle di MatLab. La sua importanza è tale che molte altre librerie popolari, come Pandas, SciPy, scikit-learn, TensorFlow e PyTorch, si basano su di essa per le loro operazioni numeriche.

Perché scegliere NumPy

• Prestazioni: gli array NumPy sono implementati in C e risiedono in blocchi di memoria contigui, rendendo le operazioni molto più veloci rispetto alle liste Python native.
• Compatibilità: funge da base per l'intero ecosistema scientifico Python, garantendo interoperabilità tra librerie.
• Espressività: la sintassi è vicina al linguaggio matematico, consentendo di scrivere codice chiaro e compatto.

Installazione e import

Per installare NumPy è sufficiente eseguire:

pip install numpy

E per importarlo nei propri script si utilizza convenzionalmente l'alias np:

import numpy as np

L'oggetto ndarray

Il cuore di NumPy è l'array N-dimensionale, o ndarray. Si tratta di una griglia di valori, tutti dello stesso tipo (dtype), con dimensioni fisse. A differenza delle liste Python, un ndarray non supporta operazioni dinamiche come append o delete in-place: per modificare la dimensione è necessario creare un nuovo array.

Creazione di un ndarray

Esistono diversi modi per creare un array NumPy. I più comuni sono i seguenti.

Da una lista Python:

a = np.array([1, 2, 3])

Array multidimensionale:

b = np.array([[1, 2], [3, 4]])

Array di zeri e di uni:

zeros = np.zeros((3, 3))
ones = np.ones((2, 2))

Matrice identità:

identity = np.eye(3)

Sequenze lineari:

seq = np.arange(0, 10, 2)
lin = np.linspace(0, 1, 5)

seq produce l'array [0 2 4 6 8], mentre lin restituisce 5 valori equidistanti tra 0 e 1.

Proprietà fondamentali

Un ndarray espone diverse proprietà utili per conoscerne la struttura:

• ndim: numero di dimensioni.
• shape: tupla che descrive la dimensione di ogni asse.
• size: numero totale di elementi.
• dtype: tipo di dato degli elementi.
• itemsize: numero di byte occupati da ogni elemento.

Esempio su un array bidimensionale:

a2D = np.array([[1, 2, 3], [4, 5, 6]])
print('ndim:', a2D.ndim)
print('shape:', a2D.shape)
print('size:', a2D.size)
print('dtype:', a2D.dtype)
print('itemsize:', a2D.itemsize)

L'output restituisce ndim = 2, shape = (2, 3), size = 6, dtype = int64 e itemsize = 8 byte per elemento sui sistemi 64-bit più diffusi come Linux e macOS. Su Windows il tipo intero di default può invece risultare int32: è quindi buona prassi non dare per scontato il dtype predefinito e, quando la portabilità del codice tra sistemi operativi diversi è importante, specificarlo sempre esplicitamente in fase di creazione dell'array.

Per un array tridimensionale:

a3D = np.array([[[1, 2, 3], [4, 5, 6]], [[7, 8, 9], [10, 11, 12]]])
print('ndim:', a3D.ndim)
print('shape:', a3D.shape)
print('size:', a3D.size)

In questo caso ndim vale 3, shape è (2, 2, 3) e size è 12.

Shape e reshape, la flessibilità della vista

Gli array NumPy sono allocati in un unico blocco di memoria contiguo, anche quando sono multidimensionali. Questo garantisce accessi veloci e permette di separare il data buffer (i dati grezzi) dalla view (l'interpretazione di quei dati). La forma (shape) di un array può essere modificata senza alterare il buffer sottostante, semplicemente creando una nuova vista.

Partiamo da un array monodimensionale di 12 elementi:

a = np.arange(0, 12)
print(a.shape)

Possiamo ridefinirlo come matrice 3x4:

b = a.reshape((3, 4))
print(b.shape)
print(b)

Il risultato è la matrice seguente.

[[ 0  1  2  3]
 [ 4  5  6  7]
 [ 8  9 10 11]]

Oppure come vettore colonna (12x1):

c = a.reshape((12, 1))
print(c.shape)

In entrambi i casi b e c condividono lo stesso buffer di a. Questa condivisione è alla base del concetto di vista, che approfondiamo nella prossima sezione.

Layout di memoria, ordine C e ordine Fortran

Il modo in cui gli elementi di un array multidimensionale vengono effettivamente disposti nel blocco di memoria contiguo dipende da un parametro chiamato ordine di memorizzazione. NumPy supporta due convenzioni principali:

• Ordine C (row-major): gli elementi di ogni riga sono contigui in memoria. È l'ordine di default in NumPy e coincide con la convenzione del linguaggio C.
• Ordine Fortran (column-major): gli elementi di ogni colonna sono contigui in memoria. È la convenzione usata da Fortran e da alcune librerie di algebra lineare storiche.

Questo diventa rilevante quando NumPy scambia dati con librerie esterne che si aspettano un ordine specifico, oppure quando si vuole ottimizzare l'accesso sequenziale a righe o colonne per ridurre i cache miss.

a = np.array([[1, 2, 3], [4, 5, 6]], order='C')
b = np.array([[1, 2, 3], [4, 5, 6]], order='F')

Anche reshape() accetta lo stesso parametro, per controllare l'ordine con cui gli elementi vengono riletti durante la ridefinizione della forma:

c = np.arange(6).reshape((2, 3), order='F')

Per il lettore alle prime armi la regola pratica è: lasciare il default (order='C') a meno che non si stia interfacciando NumPy con codice esterno che richiede esplicitamente l'ordine Fortran.

Viste e copie, quando i dati sono condivisi

Una delle caratteristiche più importanti, e talvolta insidiose, di NumPy è che lo slicing di un ndarray crea una vista, non una copia. Ciò significa che due variabili diverse possono riferirsi allo stesso data buffer: di conseguenza, modifiche effettuate su una vista si riflettono su tutte le altre viste che condividono il buffer.

Il motivo principale per cui NumPy privilegia le viste è prestazionale: non viene allocata nuova memoria e le operazioni risultano più veloci.

a = np.arange(0, 4)
b = a[:]
b[0] = 9
print(a[0])

In questo esempio anche a[0] diventa 9, perché b è una vista sullo stesso buffer di a. Per evitare questo comportamento è necessario creare una copia esplicita tramite il metodo copy():

a = np.arange(0, 4)
b = a[:].copy()
b[0] = 9
print(a[0])

Qui a[0] resta 0, perché b è ora un array indipendente.

Slicing e indexing, un discrimine cruciale

In NumPy l'accesso tramite slicing (con la notazione dei due punti) crea sempre una vista, mentre l'accesso tramite indici interi, liste di indici (fancy indexing) o maschere booleane crea sempre una copia.

arr = np.array([1, 2, 4, 8, 16, 32])
vista = arr[1:4:2]
print(vista.base)
print(vista.flags.owndata)
copia = arr[[1, 3]]
print(copia.base)
print(copia.flags.owndata)
mask = [False, True, False, True, False, False]
copia_bool = arr[mask]

La tabella seguente riassume il comportamento delle diverse tecniche di accesso.

Questo comportamento ha implicazioni pratiche: se si modifica l'array originale dopo aver creato una vista, la vista riflette il cambiamento, mentre la copia rimane invariata.

arr[1] = 100
print(vista)
print(copia)

Dopo la modifica, vista diventa [100 8], mentre copia resta [2 8].

Universal Functions, operazioni vettorizzate

Le operazioni elemento per elemento (element-wise) sono all'ordine del giorno in data science. Con le liste Python, per eseguirle è necessario ciclare esplicitamente con dei loop, il che risulta estremamente lento per grandi volumi di dati.

NumPy risolve il problema con le Universal Functions (ufunc): funzioni che operano su tutti gli elementi di un array in modo vettorizzato, eseguendo il calcolo direttamente in C. NumPy mette a disposizione numerose decine di ufunc, e molte di esse sono accessibili anche tramite overload degli operatori, rendendo il codice più leggibile.

x1 = np.array([1, 2, 3])
x2 = np.array([4, 5, 6])
somma = np.add(x1, x2)
differenza = np.subtract(x1, x2)
prodotto = np.multiply(x1, x2)
divisione = np.divide(x1, x2)
potenza = np.power(x1, 2)
radice = np.sqrt(x1)
matrice_prodotto = np.matmul(x1.reshape(1, 3), x2.reshape(3, 1))

Le prime cinque operazioni sono equivalenti, rispettivamente, agli operatori +, -, *, / e **. L'ultima riga, np.matmul, esegue invece un prodotto matriciale vero e proprio, e non un'operazione element-wise.

Un test di prestazioni

Per apprezzare la differenza di velocità, si consideri l'aggiunta di 1 a ogni elemento di un array di 10 milioni di elementi.

import numpy as np
n = 10000000
arr = np.arange(n, dtype=np.int64)
lst = list(range(n))
t_numpy = arr + 1
t_list_comp = [x + 1 for x in lst]
t_list_loop = lst[:]
for i in range(len(t_list_loop)):
    t_list_loop[i] += 1

I tempi medi rilevati su un sistema di riferimento sono i seguenti:

• NumPy: 0,005304 s
• List comprehension: 0,303015 s
• Loop esplicito: 0,553153 s

NumPy risulta essere circa due ordini di grandezza più veloce, grazie alla vettorizzazione e alla gestione contigua della memoria.

Broadcast, quando le forme non corrispondono

Cosa succede quando si applica una ufunc tra due array di forma diversa, o tra un array e uno scalare? NumPy cerca di riadattare gli array più piccoli attraverso un meccanismo chiamato broadcast, che espande virtualmente le dimensioni mancanti o di dimensione 1 per renderli compatibili.

Regole del broadcast

1. Allineamento delle dimensioni: se i due array hanno un numero diverso di dimensioni, quello con meno dimensioni viene allungato a sinistra aggiungendo assi di dimensione 1. Esempio: a.shape = (2, 3, 4), b.shape = (3, 4) diventa b trattato come (1, 3, 4).
2. Confronto dimensione per dimensione, partendo da destra: se le dimensioni sono uguali l'operazione è valida; se una delle due è 1, quell'asse viene allungato per eguagliare l'altro senza copiare fisicamente i dati; se entrambe sono diverse e nessuna è 1, si genera un errore.

Esempi di broadcast

a = np.array([1.0, 2.0, 3.0])
b = 2.0
c = a * b
print(c)

In questo caso b viene trattato come [2.0, 2.0, 2.0] e il risultato è [2.0, 4.0, 6.0].

a = np.arange(12).reshape(3, 4)
b = np.array([10, 20, 30, 40])
c = a + b

Qui b, di forma (4,), viene allungato a (3, 4) sommando riga per riga.

Casi validi di broadcast:

• a.shape = (2, 3, 4), b.shape = (1, 3, 4): valido, b viene allungato sul primo asse.
• a.shape = (2, 1, 4), b.shape = (1, 3, 4): valido, entrambi vengono allungati per ottenere (2, 3, 4).

Caso di errore:

• a.shape = (2, 2, 4), b.shape = (1, 3, 4): errore, perché la seconda dimensione (2 contro 3) non è compatibile.

Il risultato di un'operazione con broadcast ha forma pari alla massima tra le dimensioni di input per ogni asse.

Creazione avanzata di array

Oltre alla conversione da liste, NumPy offre routine dedicate per creare array con valori significativi e per comporre array esistenti.

Array creation routines

x = np.zeros((3, 2))
x = np.ones(4)
x = np.empty((4, 3))

Va precisato che np.empty() non genera valori casuali: restituisce un array la cui memoria non viene inizializzata, contenente quindi qualunque dato fosse già presente in quella porzione di RAM al momento dell'allocazione. È più veloce di np.zeros() proprio perché salta la fase di azzeramento, ma va usato solo quando si ha la certezza di sovrascrivere ogni elemento prima di leggerlo.

Composizione di array esistenti

È possibile impilare array verticalmente (vstack) o orizzontalmente (hstack).

Stack verticale, per righe:

a1 = np.array([1, 2, 3])
a2 = np.array([4, 5, 6])
v = np.vstack((a1, a2))

Il risultato è la matrice [[1 2 3] [4 5 6]].

Stack orizzontale, per colonne:

h = np.hstack((a1, a2))

Il risultato è l'array [1 2 3 4 5 6].

Modifica di un ndarray

Le operazioni di modifica, come inserimento, cancellazione e concatenazione, non avvengono in-place ma restituiscono un nuovo array. È importante ricordarlo per evitare sprechi di memoria o comportamenti inattesi.

Inserimento con np.insert

a = np.array([[1, 1], [2, 2], [3, 3]])
np.insert(a, 1, 5, axis=1)

Il risultato inserisce una colonna di valori 5 tra le colonne esistenti, ottenendo [[1 5 1] [2 5 2] [3 5 3]].

Concatenazione con np.concatenate

Concatenazione per righe (axis=0):

a = np.array([[1, 2], [3, 4]])
b = np.array([[5, 6]])
np.concatenate((a, b), axis=0)

Concatenazione per colonne (axis=1), che richiede b trasposto:

np.concatenate((a, b.T), axis=1)

Concatenazione su un unico asse, con appiattimento del risultato:

np.concatenate((a, b), axis=None)

Append e delete

arr = np.array([[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12]])
np.delete(arr, 1, 0)

Questa istruzione rimuove la seconda riga (indice 1) lungo l'asse 0, restituendo [[1 2 3 4] [9 10 11 12]].

Strumenti avanzati: einsum e Numba

Per chi supera le operazioni di base, NumPy e il suo ecosistema offrono due strumenti particolarmente utili quando le ufunc standard non bastano più.

np.einsum per operazioni tensoriali

np.einsum permette di esprimere operazioni tra array usando la notazione degli indici tipica del calcolo tensoriale, offrendo un controllo fine su quali assi vengono sommati o mantenuti. È particolarmente utile quando si lavora con tensori multidimensionali e si vuole evitare la creazione di array intermedi inutili, con benefici sia in termini di leggibilità del codice sia, in molti casi, di prestazioni.

a = np.arange(6).reshape(2, 3)
b = np.arange(6).reshape(3, 2)
prodotto = np.einsum('ij,jk->ik', a, b)

La stringa 'ij,jk->ik' descrive esplicitamente quali indici vengono moltiplicati e sommati, rendendo np.einsum uno strumento potente ma da usare con attenzione, poiché una notazione errata produce risultati silenziosamente sbagliati anziché un errore.

Numba per l'accelerazione dei loop

Quando un algoritmo non si presta facilmente alla vettorizzazione, ad esempio per la presenza di dipendenze sequenziali tra un'iterazione e la successiva, Numba permette di compilare funzioni Python in codice macchina al primo utilizzo (compilazione Just-In-Time), avvicinando le prestazioni a quelle del C senza abbandonare la sintassi Python.

from numba import jit
@jit(nopython=True)

def somma_cumulativa(arr):

risultato = np.zeros_like(arr)

acc = 0

for i in range(len(arr)):

acc += arr[i]

risultato[i] = acc

return risultato

L'opzione nopython=True forza Numba a compilare la funzione senza ricadere sull'interprete Python standard in caso di codice non supportato, garantendo il massimo guadagno prestazionale.

Best practice per un uso efficiente di NumPy

Per ottenere il massimo dalle prestazioni e dalla chiarezza del codice, si consiglia di:

• Specificare sempre il dtype quando si creano array, per controllare l'occupazione di memoria (ad esempio np.float32 invece di float64 se la precisione non è critica) e per garantire un comportamento coerente tra sistemi operativi diversi.
• Preferire np.asarray() a np.array() quando si ricevono dati che potrebbero già essere un ndarray: np.asarray() evita di creare una copia se l'input è già un array NumPy.
• Evitare loop Python espliciti: utilizzare sempre le ufunc vettorizzate o, se strettamente necessario, strumenti come Numba per compilare il loop.
• Utilizzare il broadcasting invece di espandere manualmente gli array con np.tile o np.repeat, per risparmiare memoria e tempo.
• Sfruttare l'indicizzazione avanzata (fancy e boolean) per selezioni complesse, ricordando però che producono sempre copie e non viste.
• Per operazioni tensoriali complesse, valutare l'uso di np.einsum, che offre un controllo fine sulla contrazione degli indici.
• Prestare attenzione all'ordine di memorizzazione (order='C' o order='F') quando si interfaccia NumPy con librerie esterne che si aspettano una convenzione specifica.

NumPy è uno strumento potente e flessibile. Comprendere a fondo i meccanismi di memoria, vista e broadcast è fondamentale per scrivere codice performante e privo di effetti collaterali. Con queste basi, si è pronti ad affrontare con sicurezza l'intero ecosistema scientifico Python.