6.13. Algebră liniară

Algebra liniară pe o cameră înseamnă lucrul cu matrice mici: o rotație 3x3 care fuzionează un eșantion IMU în cadrul lumii, o matrice de calibrare care corectează o lentilă, actualizarea de stare-covarianță a unui filtru Kalman, o potrivire polinomială ale cărei ecuații normale se reduc la o rezolvare liniară minusculă. Submodulele numpy.linalg și scipy.linalg acoperă exact această scară.

Funcțiile se află în două module. Operațiile de produs matricial se găsesc la nivelul superior numpy; descompunerile și inversa matricei sunt în numpy.linalg; rezolvitoarele dedicate de sisteme liniare sunt în scipy.linalg. O înmulțire de matrice 2 pe 2, de exemplu:

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. Ce este disponibil

  • dot() – produs matricial sau vectorial.

  • cross() – produs vectorial 3-D.

  • trace() – suma diagonalei.

  • inv() – inversa unei matrice.

  • det() – determinant.

  • cholesky() – descompunere Cholesky (intrare simetrică pozitiv definită).

  • eig() – valori proprii și vectori proprii ai unei matrice reale simetrice.

  • norm() – norma 2 a unui vector sau a unei matrice.

  • qr() – descompunere QR cu mode='reduced' (implicit) sau mode='complete'.

  • solve_triangular() – rezolvă A @ x = b când A este triunghiulară.

  • cho_solve() – rezolvă A @ x = b dat fiind un factor Cholesky al lui A.

6.13.2. dot, cross, trace

dot() este modul în care se scrie înmulțirea de matrice pe cameră. Operatorul @ pe care numpy de pe desktop îl oferă pentru aceeași sarcină nu este implementat pe ndarray, așa că fiecare produs matrice-vector, matrice-matrice și produs scalar de vectori trece prin apelul 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

Rezultatul lui dot() este întotdeauna de dtype float. cross() este produsul vectorial a doi vectori 3-D, iar trace() suma diagonalei principale a unei matrice pătratice.

6.13.3. inv și 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))

Inversa este calculată prin eliminare Gauss-Jordan, așadar inv() ridică ValueError atunci când matricea este singulară (o intrare diagonală devine zero în timpul eliminării). Costul de RAM este aproximativ dublul dimensiunii intrării.

Determinantul reutilizează aceeași eliminare – timpul de execuție este în esență același cu cel al inversei.

Când scopul aplicației este rezolvarea unui sistem liniar, nu inversa și înmulți – preferă rezolvitoarele dedicate de mai jos. Ambele sunt mai rapide și se comportă mai bine din punct de vedere numeric.

6.13.4. cholesky

Pentru o matrice simetrică pozitiv definită A, cholesky() returnează un factor inferior triunghiular L astfel încât A = L @ L.T

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

Dacă intrarea nu este pozitiv definită sau nu este simetrică, se ridică ValueError.

Factorul Cholesky reprezintă jumătate din efortul unei factorizări LU și este punctul de plecare potrivit pentru orice problemă a cărei matrice se știe că este simetrică pozitiv definită (actualizări de covarianță, ecuații normale dintr-o potrivire prin metoda celor mai mici pătrate).

6.13.5. eig

eig() funcționează doar pe matrice reale simetrice. Matricele nesimetrice ridică ValueError. Returnează un tuplu de 2 elemente (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:

  • Valorile proprii sunt returnate fără o ordine anume. Aplică sort() (și aceeași permutare asupra vectorilor proprii prin argsort()) când ordinea sortată contează.

  • Un vector propriu este unic doar până la un scalar nenul, așa că semnul vectorilor proprii individuali nu este definit în mod unic. Două rulări corecte pot produce vectori de semn opus; acest lucru este inofensiv.

6.13.6. norm

Norma euclidiană (Frobenius) a unui vector sau a unei 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...

Cuvântul-cheie opțional axis= calculează norma de-a lungul unei singure axe, în loc de pe întregul tablou.

6.13.7. qr

qr() descompune o matrice dreptunghiulară A (formă (M, N)) într-o matrice ortonormală Q și o matrice superior triunghiulară R astfel încât 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)

Descompunerea este implementată prin rotații Givens succesive. Alegerea potrivită pentru problemele de cele mai mici pătrate în care matricea nu este simetrică.

6.13.8. Rezolvarea sistemelor

Cele două rezolvitoare dedicate din ulab.scipy.linalg sunt deopotrivă mai rapide și mai precise decât np.dot(np.linalg.inv(A), b):

  • solve_triangular(a, b, lower=False)() – rezolvă a @ x = b presupunând că a este triunghiulară:

    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)() – dat fiind un factor Cholesky L, rezolvă A @ x = b unde A = L @ L.T

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

Recurge la acestea în loc de a inversa ori de câte ori structura lui A o permite – ele economisesc efortul de eliminare și inversa explicită.

6.13.9. Rezolvarea unui sistem liniar mic

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

Aceeași problemă poate fi exprimată prin factorizarea lui A și apelarea rezolvitorului potrivit – mai rapid și mai precis atunci când A are structura corespunzătoare.

Pentru referința completă la nivel de argumente, vezi numpy.linalg — Rutine de algebră liniară și scipy.linalg — Rutine de algebră liniară.