Python per l’algebra lineare: numpy

Python per l’algebra lineare: numpy#

La Numerical Linear Algebra (NLA) è la materia che si occupa di studiare come le operazioni tra matrici e vettori possono essere utilizzate all’interno di algoritmi con lo scopo di produrre risposte efficienti e accurate a diversi problemi: da quelli di carattere matematico (come la risoluzione approssimata di equazioni differenziali), a quelli di carattere ingegneristico (come la distribuzione ottimale delle risorse, o le predizioni di outcome futuri dai dati).

Di conseguenza, è di vitale importanza avere un modo efficiente per implementare operazioni matriciali e vettori. Nel seguito (come è comune fare anche nel codice), indicheremo con le lettere maiuscole (A, B, …) le matrici, con lettere minuscole (v, w, …) i vettori.

Come abbiamo già accennato a lezione, la libreria Python più comunemente utilizzata per lavorare con matrici è vettori è numpy.

Durante il corso, faremo largo uso di numpy, nonché di alcune librerie add-on (come scipy e matplotlib).

Nel caso in cui numpy non sia ancora installato sul vostro virtual environment, si ricorda di eseguire i seguenti comandi per installarlo:

conda activate CN24
conda install numpy

Note

Durante il corso utilizzeremo varie funzioni di numpy, ma chiaramente ne esistono molte altre. Una lista completa di tutte le funzioni, oltre che una documentazione sul loro utilizzo, può essere trovata sul sito ufficiale numpy.org.

Creazione di array con numpy#

numpy implementa un tipo di dato speciale che estende i classici array di Python (come liste, tuple, …), specificatamente pensato per massimizzare l’efficienza di operazioni di algebra lineare, chiamati ndarray. Tramite gli ndarray, è possibile inizializzare facilmente vettori, matrici, tensori, …

Questi non solo avranno tutte le funzionalità degli altri array di Python, ma implementeranno in maniera nativa tutte le operazioni tipiche tra vettori e matrici, come la somma elemento per elemento, il prodotto riga-per-colonna, e simili.

Unico limite degli ndarray è che essi possono contere, al loro interno, solo numeri (e non stringhe o altro, come per le liste o le tuple).

Il metodo classico per la creazione di un ndarray è quello di castarlo da un’array classico di Python. Vediamo come creare, ad esempio, un vettore:

import numpy as np

# Creo un vettore numpy
a: list[int] = [1, 2, 3]
a_vec: np.ndarray = np.array(a)
print(type(a_vec)) # Verifichiamo il tipo di dato di a_vec
print(a_vec)

<class 'numpy.ndarray'>
[1 2 3]

Una delle proprietà di base degli ndarray è la propria shape, che rappresenta la loro dimensione. Più precisamente, un vettore avrà come shape una tupla contenente un singolo elemento, che ne rappresenta la lunghezza, una matrice avrà come shape una tupla contenente due elementi: il numero di righe e il numero di colonne, rispettivamente.

Per inizializzare una matrice a partire da un’array Python, è necessario costruire un’array composto da altri array, ciascuno dei quali rappresenta una riga della matrice che si vuole creare.

Warning

Per ovvie ragioni, è necessario che l’array che definisce la matrice abbia tutte le righe della stessa lunghezza. In caso contrario, il codice restituirà errore.

Vediamo ad esempio come costruire una matrice e un vettore, e stamparne la shape.

# Definiamo un vettore v e una matrice A
v = np.array([1, 1, 3, 2])
A = np.array([[1, 1, -1],
              [0, -1, 1]])

# Stampiamo la shape
print(f"La shape di v è: {v.shape}.")
print(f"La shape di A è: {A.shape}.")

La shape di v è: (4,).
La shape di A è: (2, 3).

Altra proprietà chiave degli ndarray è il loro dtype. Infatti, proprio per il fatto che sono stati costruiti per eseguire operazioni di algebra lineare, numpy implementa vari tipi differenti di aritmetica floating point per i propri array. E’ infatti possibile definire array a precisione singola (np.float32), doppia (np.float64) o anche mezza (np.float16). Questa informazione è rappresentata, appunto, dal dtype.

Se non specificato in fase di creazione, gli ndarray hanno dtype = np.float64 se contengono numeri floating point, np.int64 altrimenti.

# Creiamo un vettore
v = np.array([1, 2, 1.2])
print(v.dtype) # np.float64

# Creiamo lo stesso vettore in precisione singola
v = np.array([1, 2, 1.2], dtype=np.float32)
print(v.dtype)

float64
float32

Funzioni per la creazione di numpy array#

Nella maggior parte delle applicazioni pratiche, le matrici e i vettori utilizzati sono di dimensioni enormi, dell’ordine di migliaia se non di milioni di elementi. Ovviamente, questo rende impossibile definirli a partire da array di Python, come visto precedentemente.

Fortunatamente, quando gli array che si vuol generare possiedono un pattern specifico, è possibile inizializzarle tramite una serie di funzioni, riportate sotto.

  • np.linspace(a, b, n): Crea un vettore di lunghezza n, che contiene n elementi uniformemente distributi nel’intervallo \([a, b]\), estremi inclusi.

  • np.arange(inizio, fine, passo): Simile alla funzione Python range. Crea un vettore che contiene tutti i numeri interi a partire da inizio, fino a fine-1, con passo fissato.

  • np.zeros((m, n)): Crea una matrice di dimensione (m, n) di zeri. Chiaramente, per creare un vettore invece che una matrice, si utilizza la funzione np.zeros((m, )).

  • np.ones((m, n)): Come prima, ma crea una matrice (o vettore) di 1.

  • np.zeros_like(a): Crea una matrice o un vettore, della stessa dimensione di un altro array a. E’ equivalente a np.zeros(a.shape). Esiste l’equivalente np.ones_like(a).

  • np.diag(v): Dato un vettore v di lunghezza n, costruisce una matrice diagonale di dimensione (n, n), la cui diagonale è v.

  • np.random.randn(m, n): Crea una matrice di dimensione (m, n) di numeri casuali distribuiti con distribuzione normale standard.

  • np.random.rand(m, n): Uguale a prima, ma con valori estratti da una distribuzione uniforme nel dominio \([0, 1]\).

Costruiamo, ad esempio, un vettore di valori casuali in \([0, 1]\) di lunghezza 10.

# Creiamo il vettore richiesto
v = np.random.rand(10)
print(v)

[0.53860875 0.13699224 0.13390741 0.88724061 0.19142975 0.60958048
 0.06971019 0.30655518 0.78817803 0.77595642]

Operazioni tra ndarray#

Una volta imparato come creare degli array, possiamo passare a come utilizzarli. Infatti, come accennato precedentemente, gli ndarray implementano nativamente le operazioni classiche di algebra lineare. In particolare, quasi tutte le operazioni standard sono applicate elemento per elemento:

# Creiamo due ndarray
a = np.array([1, -1, 0])
b = np.linspace(1, 3, 3) # array(1, 2, 3)

# Eseguiamo operazioni tra loro
s = a + b
d = a - b
p = a * b
f = a / b

print(f"a = {a}, b = {b} \nSomma: {s} \nDifferenza: {d} \nProdotto: {p} \nDivisione: {f}.")

a = [ 1 -1  0], b = [1. 2. 3.] 
Somma: [2. 1. 3.] 
Differenza: [ 0. -3. -3.] 
Prodotto: [ 1. -2.  0.] 
Divisione: [ 1.  -0.5  0. ].

Similmente, lavorano elemento per elemento anche tutte le funzioni base di matematica, come l’operazione di seno, coseno, tangente, esponenziale, logaritmo, ecc.

# Definiamo un array
x = np.linspace(0, 2*np.pi, 4)

# Calcoliamone i valori trigonometrici
print(f"sin(x) = {np.sin(x)}.")
print(f"cos(x) = {np.cos(x)}.")
print(f"tan(x) = {np.tan(x)}.")

# Ed esponenziale e logaritmo
print(f"e^x = {np.exp(x)}.")
print(f"ln(x) = {np.log(x + 1)}.")
print(f"log_10(x) = {np.log10(x + 1)}.")

sin(x) = [ 0.00000000e+00  8.66025404e-01 -8.66025404e-01 -2.44929360e-16].
cos(x) = [ 1.  -0.5 -0.5  1. ].
tan(x) = [ 0.00000000e+00 -1.73205081e+00  1.73205081e+00 -2.44929360e-16].
e^x = [  1.           8.1205274   65.9429652  535.49165552].
ln(x) = [0.         1.12959244 1.64650057 1.98556831].
log_10(x) = [0.         0.49057577 0.71506611 0.86232136].

Unica, importante differenza alla regola delle operazioni elemento per elemento è rappresentato dal prodotto riga-per-colonna tra matrici e vettori (anche detto prodotto scalare o prodotto interno quando applicato a due vettori).

Si ricorda che, presa una matrice \(A = (a_{i, j})\) di dimensione \(m \times n\), e un vettore \(x = (x_1, \dots, x_n) \in \mathbb{R}^n\), il prodotto riga-per-colonna tra \(A\) e \(x\) è definito come:

\begin{align} \left( Ax \right)i = \sum{j = 1}^n a_{i, j} x_j. \end{align}

Questa operazione è implementata su numpy tramite l’operatore @ o, equivalentemente, tramite la funzione np.dot(A, x).

# Definiamo una matrice A e un vettore x
A = np.array([[1, 1, 1], [0, -1, 0], [0, 0, 1]])
x = np.array([1, 0, 1])

# Calcoliamo y = Ax
y = A@x
print(y)

[2 0 1]

E per quanto riguarda i vettori? E’ importante fare una distinzione tra tre tipologie differenti di vettori rappresentabili come ndarray. Infatti, un vettore può assumere diverse forme:

  1. un vettore classico, a dimensione singola, la cui shape è qualcosa del tipo (n,).

  2. un vettore colonna, in cui si vede il vettore come fosse una matrice con una sola colonna, e un numero di righe pari alla sua lunghezza. In questo caso, il vettore avrà una shape della forma (n, 1).

  3. un vettore riga, in cui si vede il vettore come fosse una matrice con una sola riga, e un numero di colonne pari alla sua lunghezza. In questo caso, il vettore avrà una shape della forma (1, n).

E’ importante stare attenti a questa distinzione, poiché alcune funzionalità di numpy lavorano diversamente se il vettore è un vettore classico, oppure un vettore colonna. Ad esempio, l’operatore @ restituirà il prodotto scalare standard tra vettori nel caso in cui i due elementi moltiplicati siano due vettori classici, restituirà errore nel caso in cui siano due vettori colonna, e darà il prodotto scalare standard se i vettori saranno il primo un vettore riga, il secondo un vettore colonna.

# Definiamo due vettori classici
v = np.array([1, -1, 1])
w = np.array([0, 1, 1])
print(v.shape, w.shape) # Controlliamo che siano vettori classici

# Calcoliamone il prodotto scalare
print(v @ w)

# Definiamo la loro versione come vettori colonna
v = np.array([[1], [-1], [1]])
w = np.array([[0], [1], [1]])
print(v.shape, w.shape)

# Calcoliamone il prodotto scalare
print(v @ w)

(3,) (3,)
0
(3, 1) (3, 1)

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[8], line 15
     12 print(v.shape, w.shape)
     14 # Calcoliamone il prodotto scalare
---> 15 print(v @ w)

ValueError: matmul: Input operand 1 has a mismatch in its core dimension 0, with gufunc signature (n?,k),(k,m?)->(n?,m?) (size 3 is different from 1)

Operazioni logiche tra vettori#

Sfruttando la caratteristica di Python per cui le variabili booleane sono sostanzialmente numeri (0 per il False, 1 per il True), è possibile anche definire degli ndarray che contengono valori booleani. Questi saranno particolarmente utili specialmente per svolgere in maniera efficiente operazioni di filtraggio di elementi.

Sarà possibile chiaramente eseguire operazioni booleane classiche (come and, or, not) e di confronto (==, !=, ecc.) tra ndarray booleani, che verranno applicate elemento per elemento. Di seguito una tabella che rappresenta le principali operazioni possibili tra ndarray booleani.

Operatore

Significato

==

UGUAGLIANZA

!=

DISUGUAGLIANZA

>, >=

MAGGIORE

<, <=

MINORE

&

AND

|

OR

!

NOT

 
# Definiamo tre vettori casuali
a = np.random.rand(10)
b = np.random.rand(10)
c = np.random.rand(10)

# Generiamo il vettore booleano `v` che vale True quando un elemento di `a`
# è maggiore o uguale del corrispettivo elemento di `b`
v = a >= b
print(v)

# E il vettore `w`che vale True quando un elemento di `b` è maggiore o uguale
# del corrispettivo elemento di `c`
w = b >= c
print(w)

# Ora uniamoli con un operazione di `and` elemento per elemento
print(w & v)

[False False False False False False  True  True  True False]
[ True  True  True  True False  True  True False False False]
[False False False False False False  True False False False]

Slicing su ndarray#

Chiaramente, anche sugli ndarray è possibile definire le operazioni di slicing, che verranno come di consueto effettuate tramite la parentesi quadra []. Sui vettori classici, lo slicing funziona esattamente come per liste e tuple, in cui la notazione v[inizio:fine] restituisce gli elementi di v da inizio a fine-1.

# Definiamo un array
v = np.array([0, 1, -1, 2, 1, -1])

# Slicing
w = v[0:3]
print(w)

[ 0  1 -1]

Il vero vantaggio degli ndarray rispetto alle liste e tuple in termini di slicing, è la possibilità di utilizzare gli ndarray booleani per filtrare gli elementi di un altro ndarray. Infatti, se all’interno delle parentesi quadre inseriamo un array booleano della stessa dimensione dell’array su cui stiamo facendo slicing, questo ci ritornerà tutti gli elementi dell’array che corrispondono alle posizioni dei True.

# Creiamo un array di esempio
v = np.array([0, 1, -1, 2, 1, -1])
print(v)

# Definiamo un vettore booleano
b = np.array([True, True, False, True, False, False])
print(b)

# Slicing
print(v[b])

[ 0  1 -1  2  1 -1]
[ True  True False  True False False]
[0 1 2]

Questo ci permette di eseguire in maniera estramamente rapida e in un linguaggio molto vicino a quello naturale, operazioni condizionali su vettori.

Vediamo ad esempio come prendere un vettore casuale, e porre a 0 tutti i suoi valori negativi (operazione, nota come proiezione sull’asse positivo, molto comune in alcuni problemi di algebra lineare).

# Definiamo un vettore casuale
x = np.random.randn(8)
print(x)

# Proiettiamo sull'asse positivo
x[x < 0] = 0
print(x)

[-0.80166815  0.93414683  0.62817759 -0.21499552 -0.22965556 -0.5291032
 -0.85889971 -0.30573981]
[0.         0.93414683 0.62817759 0.         0.         0.
 0.         0.        ]

Slicing di matrici#

Sulle matrici (come anche sugli ndarray di dimensione più alta) lo slicing funziona in maniera naturale, con la sola differenza che all’interno della parentesi quadra sarà necessario inserire due indici (o comandi della forma inizio:fine:step): uno per gli indici di riga, e uno per gli indici di colonna.

Vediamo ad esempio come estrarre, da una matrice A di shape (3, 3), la sua sottomatrice principale di ordine 2, ovvero la sottomatrice di dimensione (2, 2) che si trova nell’angolo in alto a sinistra di A.

# Creiamo la matrice
A = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])

# Slicing
B = A[:2, :2]
print(B)

[[1 2]
 [4 5]]

Manipolazione di matrici e vettori#

Ultima proprietà che vediamo degli ndarray, sono tutte le operazioni comunemente effettuate nell’algebra lineare su di essi, come ad esempio il calcolo della norma, dell’inversa, del rango, della trasposta, oltre che di operazioni di utilità come la funzione reshape, che modifica la forma di un vettore o matrice, senza cambiarne il contenuto.

Sotto un elenco di tali operazioni. Notare come, alcune di esse, richiedono di accedere ad un sottopacchetto di numpy, ovvero np.linalg, che come suggerisce il nome, implementa varie operazioni di algebra lineare.

  • np.linalg.norm(a, p): Calcola la norma-\(p\) di un vettore o di una matrice.

  • np.linalg.cond(A, p): Calcola il numero di condizionamento di una matrice in norma \(p\);

  • np.linalg.matrix_rank(A): Calcola il rango della matrice A;

  • np.linalg.inv(A): Calcola l’inverso della matrice A, quando esiste. ATTENZIONE: Questa operazione è molto lenta anche per matrici relativamente piccole;

  • np.transpose(A): Calcola la trasposta della matrice A. E’ equivalente a A.T;

  • np.reshape(a, new_shape): Modifica la shape di un ndarray a, nella forma specificata.

Vediamo un esempio di applicazione di tali funzioni.

# Definiamo un vettore v
v = np.linspace(1, 9, 9) # (1, 2, 3, ..., 9)
print(v)

# E costruiamo la matrice A di shape 3x3 ottenuta modificando v
A = np.reshape(v, (3, 3))
print(A)

# Controlliamo se A è di rango massimo
print(np.linalg.matrix_rank(A))

# E, se ha rango massimo, calcoliamone l'inversa
if np.linalg.matrix_rank(A) == A.shape[0]:
    A_inv = np.linalg.inv(A)

[1. 2. 3. 4. 5. 6. 7. 8. 9.]
[[1. 2. 3.]
 [4. 5. 6.]
 [7. 8. 9.]]
2

Approfondimento: efficienza di numpy#

L’obiettivo di questa sezione di approfondimento è verificare in modo empirico quanto numpy è più efficiente del Python classico per operazioni su vettori. Per farlo, andremo ad utilizzare la libreria time per misurare il tempo richiesto a Python base e numpy per svolgere le stesse operazioni.

In questo esempio, utilizzeremo un vettore di lunghezza 10_000_000 (10 milioni), casuale, e proveremo a calcolare la somma del valore assoluto dei suoi elementi (operazione nota come norma 1 del vettore, e solitamente indicata con \(|| v ||_1\)), utilizzando Python classico (attraverso due algoritmi differenti) e le funzioni di numpy.

import time
import numpy as np
import math

# Costruiamo il vettore
v = np.random.randn(10_000_000)

# || v ||_1 con Python modo 1
start_time = time.time()
norma_1 = 0
for i in range(len(v)):
    norma_1 = norma_1 + abs(v[i])
end_time = time.time()
print(f"Tempo impiegato con Python modo 1: {end_time - start_time}s")

# || v ||_1 con Python modo 2
start_time = time.time()
norma_1 = sum((abs(i) for i in v)) # list comprehension
end_time = time.time()
print(f"Tempo impiegato con Python modo 2: {end_time - start_time}s")

# || v ||_1 con numpy
start_time = time.time()
norma_1 = np.sum(np.abs(v))
end_time = time.time()
print(f"Tempo impiegato con Numpy: {end_time - start_time}s")

Tempo impiegato con Python modo 1: 0.9546482563018799s

Tempo impiegato con Python modo 2: 0.5896708965301514s
Tempo impiegato con Numpy: 0.01151275634765625s

Da questo esempio appare chiaro come Numpy, avendo implementato funzione pre-compilate, è estremamente più efficiente di Python standard, ed è quindi più preposto a eseguire operazioni di algebra lineare, dove l’efficienza è importantissima.

Contents