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.