6.13. Algebra lineare#

L’algebra lineare su una camera è lavoro su piccole matrici: una rotazione 3x3 che fonde un campione della IMU nel frame del mondo, una matrice di calibrazione che rettifica una lente, l’aggiornamento di stato-covarianza di un filtro di Kalman, un fit polinomiale le cui equazioni normali si riducono a una minuscola soluzione lineare. I sottomoduli numpy.linalg e scipy.linalg coprono esattamente questa scala.

Le funzioni risiedono in due moduli. Le operazioni di prodotto tra matrici si trovano al livello superiore di numpy; le decomposizioni e l’inversa di matrice sono sotto numpy.linalg; i risolutori dedicati di sistemi lineari sono sotto scipy.linalg. Una moltiplicazione di matrici 2-per-2, ad esempio:

from ulab import numpy as np

A = np.array([[1, 2], [3, 4]], dtype=np.float)
B = np.array([[5, 6], [7, 8]], dtype=np.float)
np.dot(A, B)
# array([[19.0, 22.0],
#        [43.0, 50.0]])

6.13.1. Cosa è disponibile#

  • dot() – prodotto matriciale o vettoriale.

  • cross() – prodotto vettoriale 3-D.

  • trace() – somma della diagonale.

  • inv() – inversa di matrice.

  • det() – determinante.

  • cholesky() – decomposizione di Cholesky (input simmetrico definito positivo).

  • eig() – autovalori e autovettori di una matrice reale simmetrica.

  • norm() – norma 2 di un vettore o di una matrice.

  • qr() – decomposizione QR con mode='reduced' (predefinita) o mode='complete'.

  • solve_triangular() – risolve A @ x = b quando A è triangolare.

  • cho_solve() – risolve A @ x = b dato un fattore di Cholesky di A.

6.13.2. dot, cross, trace#

dot() è il modo in cui si scrive la moltiplicazione di matrici sulla camera. L’operatore @ che il numpy desktop fornisce per lo stesso scopo non è implementato su ndarray, quindi ogni prodotto matrice-vettore, matrice-matrice e prodotto scalare tra vettori passa per la chiamata dot()

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

np.dot(a, b)             # 32.0 (scalar product)
np.cross(a, b)           # array([-3.0, 6.0, -3.0])

m = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
np.dot(m, a)             # matrix-vector product
np.dot(m, m)             # matrix-matrix product
np.trace(m)              # 1 + 5 + 9 = 15.0

Il risultato di dot() è sempre di dtype float. cross() è il prodotto vettoriale di due vettori a 3 componenti, trace() la somma della diagonale principale di una matrice quadrata.

6.13.3. inv e det#

m = np.array([[1, 2, 3, 4],
              [4, 5, 6, 4],
              [7, 9, 9, 4],
              [3, 4, 5, 6]])
print(np.linalg.inv(m))
print(np.linalg.det(m))

L’inversa viene calcolata mediante eliminazione di Gauss-Jordan, quindi inv() solleva ValueError quando la matrice è singolare (un elemento diagonale diventa zero durante l’eliminazione). Il costo in RAM è all’incirca il doppio della dimensione dell’input.

Il determinante riutilizza la stessa eliminazione – il tempo di esecuzione è essenzialmente lo stesso dell’inversa.

Quando l’obiettivo dell’applicazione è risolvere un sistema lineare, non invertire e moltiplicare – preferisci i risolutori dedicati riportati sotto. Entrambi sono più veloci e numericamente più affidabili.

6.13.4. cholesky#

Per una matrice simmetrica definita positiva A, cholesky() restituisce una matrice triangolare inferiore L tale che A = L @ L.T

a = np.array([[25, 15, -5],
              [15, 18,  0],
              [-5,  0, 11]])
L = np.linalg.cholesky(a)

Se l’input non è definito positivo o non è simmetrico, viene sollevato ValueError.

Il fattore di Cholesky richiede metà del lavoro di una fattorizzazione LU ed è il punto di partenza giusto per qualsiasi problema la cui matrice sia nota essere simmetrica definita positiva (aggiornamenti di covarianza, equazioni normali da un fit ai minimi quadrati).

6.13.5. eig#

eig() funziona solo su matrici reali simmetriche. Le matrici non simmetriche sollevano ValueError. Restituisce una 2-tupla (eigenvalues, eigenvectors)

a = np.array([[1, 2, 1, 4],
              [2, 5, 3, 5],
              [1, 3, 6, 1],
              [4, 5, 1, 7]], dtype=np.uint8)
x, y = np.linalg.eig(a)

Note:

  • Gli autovalori vengono restituiti senza un ordine particolare. Applica sort() (e la stessa permutazione agli autovettori tramite argsort()) quando l’ordinamento è importante.

  • Un autovettore è univoco solo a meno di uno scalare non nullo, quindi il segno dei singoli autovettori non è definito in modo univoco. Due esecuzioni corrette possono produrre vettori di segno opposto; questo è innocuo.

6.13.6. norm#

La norma euclidea (di Frobenius) di un vettore o di una matrice:

v = np.array([1, 2, 3, 4, 5])
m = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])

np.linalg.norm(v)        # 7.416...
np.linalg.norm(m)        # 16.881...

La parola chiave opzionale axis= calcola la norma lungo un singolo asse invece che sull’intero array.

6.13.7. qr#

qr() fattorizza una matrice rettangolare A (forma (M, N)) in una matrice ortonormale Q e una matrice triangolare superiore R tali che A == Q @ R

A = np.arange(6).reshape((3, 2))

q, r = np.linalg.qr(A)
# mode='reduced' (default): q is (3, 2), r is (2, 2)

q, r = np.linalg.qr(A, mode='complete')
# q is (3, 3), r is (3, 2)

La decomposizione è implementata tramite successive rotazioni di Givens. È la scelta giusta per problemi ai minimi quadrati in cui la matrice non è simmetrica.

6.13.8. Risolvere i sistemi#

I due risolutori dedicati sotto ulab.scipy.linalg sono entrambi più veloci e più accurati di np.dot(np.linalg.inv(A), b):

  • solve_triangular(a, b, lower=False)() – risolve a @ x = b assumendo che a sia triangolare:

    A = np.array([[3, 0, 0, 0],
              [2, 1, 0, 0],
              [1, 0, 1, 0],
              [1, 2, 1, 8]])
b = np.array([4, 2, 4, 2])
x = sp.linalg.solve_triangular(A, b, lower=True)

  • cho_solve(L, b)() – dato un fattore di Cholesky L, risolve A @ x = b dove A = L @ L.T

    L = np.linalg.cholesky(A)
x = sp.linalg.cho_solve(L, b)

Ricorri a questi invece di invertire ogni volta che la struttura di A lo consente – risparmiano lavoro di eliminazione e l’inversa esplicita.

6.13.9. Risolvere un piccolo sistema lineare#

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

x = np.dot(np.linalg.inv(A), b)
print(x)
print(np.dot(A, x))         # should equal b

Lo stesso problema può essere espresso fattorizzando A e chiamando il risolutore appropriato – più veloce e più accurato quando A ha la struttura giusta.

Per il riferimento completo a livello di argomenti, vedi numpy.linalg — Routine di algebra lineare e scipy.linalg — Routine di algebra lineare.