6.13. Lineaire algebra

Lineaire algebra op een camera is werk met kleine matrices: een 3x3-rotatie die een IMU-sample samenvoegt met het wereldframe, een kalibratiematrix die een lens rectificeert, de update van de toestandscovariantie van een Kalman-filter, een polynoomfit waarvan de normaalvergelijkingen neerkomen op een kleine lineaire oplossing. De submodules numpy.linalg en scipy.linalg dekken precies die schaal.

De functies bevinden zich in twee modules. Matrixproductbewerkingen staan op het topniveau van numpy; decomposities en de matrixinverse staan onder numpy.linalg; de specifieke oplossers voor lineaire stelsels staan onder scipy.linalg. Een 2-bij-2-matrixvermenigvuldiging bijvoorbeeld:

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. Wat beschikbaar is

  • dot() – matrix- of vectorproduct.

  • cross() – 3D-vectorkruisproduct.

  • trace() – som van de diagonaal.

  • inv() – matrixinverse.

  • det() – determinant.

  • cholesky() – Cholesky-decompositie (symmetrische positief-definiete invoer).

  • eig() – eigenwaarden en eigenvectoren van een reële symmetrische matrix.

  • norm() – 2-norm van een vector of matrix.

  • qr() – QR-decompositie met mode='reduced' (standaard) of mode='complete'.

  • solve_triangular() – los A @ x = b op wanneer A driehoekig is.

  • cho_solve() – los A @ x = b op gegeven een Cholesky-factor van A.

6.13.2. dot, cross, trace

dot() is hoe matrixvermenigvuldiging op de camera wordt geschreven. De @-operator die desktop-numpy voor dezelfde taak biedt, is niet geïmplementeerd op ndarray, dus elk matrix-vector-, matrix-matrix- en vector-puntproduct loopt via de dot()-aanroep:

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

Het resultaat van dot() is altijd van dtype float. cross() is het kruisproduct van twee 3-vectoren, trace() de som van de hoofddiagonaal van een vierkante matrix.

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

De inverse wordt berekend met Gauss-Jordan-eliminatie, dus inv() werpt ValueError op wanneer de matrix singulier is (een diagonaalelement wordt nul tijdens de eliminatie). De RAM-kosten zijn ongeveer twee keer de grootte van de invoer.

De determinant hergebruikt dezelfde eliminatie – de uitvoeringstijd is in wezen gelijk aan die van de inverse.

Wanneer het doel van de applicatie is om een lineair stelsel op te lossen, inverteer dan niet en vermenigvuldig niet – gebruik bij voorkeur de specifieke oplossers hieronder. Beide zijn sneller en numeriek beter van gedrag.

6.13.4. cholesky

Voor een symmetrische positief-definiete matrix A retourneert cholesky() een benedendriehoekige L zodat A = L @ L.T

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

Als de invoer niet positief-definiet of niet symmetrisch is, wordt ValueError opgeworpen.

De Cholesky-factor kost de helft van het werk van een LU-factorisatie en is het juiste startpunt voor elk probleem waarvan de matrix bekend symmetrisch positief-definiet is (covariantie-updates, normaalvergelijkingen uit een kleinste-kwadratenfit).

6.13.5. eig

eig() werkt alleen op reële symmetrische matrices. Niet-symmetrische matrices werpen ValueError op. De functie retourneert een 2-tuple (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)

Opmerkingen:

  • De eigenwaarden komen in willekeurige volgorde terug. Pas sort() toe (en dezelfde permutatie op de eigenvectoren via argsort()) wanneer een gesorteerde volgorde van belang is.

  • Een eigenvector is alleen uniek op een niet-nul scalair na, dus het teken van afzonderlijke eigenvectoren is niet eenduidig gedefinieerd. Twee correcte uitvoeringen kunnen vectoren met tegengesteld teken opleveren; dit is onschadelijk.

6.13.6. norm

De Euclidische (Frobenius-)norm van een vector of matrix:

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...

Het optionele sleutelwoord axis= neemt de norm langs een enkele as in plaats van over de hele array.

6.13.7. qr

qr() factoriseert een rechthoekige matrix A (vorm (M, N)) in een orthonormale Q en een bovendriehoekige R zodat 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)

De decompositie is geïmplementeerd via opeenvolgende Givens-rotaties. De juiste keuze voor kleinste-kwadratenproblemen waarbij de matrix niet symmetrisch is.

6.13.8. Stelsels oplossen

De twee specifieke oplossers onder ulab.scipy.linalg zijn beide sneller en nauwkeuriger dan np.dot(np.linalg.inv(A), b):

  • solve_triangular(a, b, lower=False)() – los a @ x = b op, aannemende dat a driehoekig is:

    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)() – gegeven een Cholesky-factor L, los A @ x = b op waarbij A = L @ L.T

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

Grijp naar deze in plaats van te inverteren wanneer de structuur van A dat toelaat – ze besparen eliminatiewerk en de expliciete inverse.

6.13.9. Een klein lineair stelsel oplossen

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

Hetzelfde probleem kan worden uitgedrukt door A te factoriseren en de juiste oplosser aan te roepen – sneller en nauwkeuriger wanneer A de juiste structuur heeft.

Voor de volledige referentie op argumentniveau, zie numpy.linalg — Routines voor lineaire algebra en scipy.linalg — Lineaire-algebraroutines.