6.13. Linjär algebra

Linjär algebra på en kamera är arbete med små matriser: en 3x3-rotation som smälter in ett IMU-sampel i världskoordinatsystemet, en kalibreringsmatris som rätar upp ett objektiv, en Kalman-filters tillstånds-kovariansuppdatering, en polynomanpassning vars normalekvationer kommer ut som ett litet linjärt lösningsproblem. Submodulerna numpy.linalg och scipy.linalg täcker exakt den skalan.

Funktionerna finns i två moduler. Matrisproduktoperationer ligger på toppnivån i numpy; dekompositioner och matrisinversen ligger under numpy.linalg; de dedikerade lösarna för linjära system ligger under scipy.linalg. En 2-gånger-2-matrismultiplikation, till exempel:

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. Vad som finns tillgängligt

  • dot() – matris- eller vektorprodukt.

  • cross() – 3D-vektors kryssprodukt.

  • trace() – summan av diagonalen.

  • inv() – matrisinvers.

  • det() – determinant.

  • cholesky() – Cholesky-dekomposition (symmetrisk positivt definit indata).

  • eig() – egenvärden och egenvektorer för en reell symmetrisk matris.

  • norm() – 2-norm för en vektor eller matris.

  • qr() – QR-dekomposition med mode='reduced' (standard) eller mode='complete'.

  • solve_triangular() – lös A @ x = b när A är triangulär.

  • cho_solve() – lös A @ x = b givet en Cholesky-faktor av A.

6.13.2. dot, cross, trace

dot() är hur matrismultiplikation skrivs på kameran. Operatorn @ som skrivbords-numpy tillhandahåller för samma uppgift är inte implementerad på ndarray, så varje matris-vektor-, matris-matris- och vektorskalärprodukt går genom anropet 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

Resultatet av dot() är alltid av dtype float. cross() är kryssprodukten av två 3-vektorer, trace() summan av huvuddiagonalen i en kvadratisk matris.

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

Inversen beräknas med Gauss-Jordan-eliminering, så inv() ger upphov till ValueError när matrisen är singulär (en diagonalpost blir noll under elimineringen). RAM-kostnaden är ungefär dubbelt så stor som indatans storlek.

Determinanten återanvänder samma eliminering – körtiden är i princip densamma som för inversen.

När applikationens mål är att lösa ett linjärt system, invertera inte och multiplicera – föredra de dedikerade lösarna nedan. Båda är snabbare och bättre numeriskt uppförda.

6.13.4. cholesky

För en symmetrisk positivt definit matris A returnerar cholesky() en nedre triangulär L så att A = L @ L.T

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

Om indatan inte är positivt definit eller inte symmetrisk ges ValueError.

Cholesky-faktorn är hälften av arbetet jämfört med en LU-faktorisering och är rätt utgångspunkt för varje problem vars matris är känd för att vara symmetrisk positivt definit (kovariansuppdateringar, normalekvationer från en minstakvadratanpassning).

6.13.5. eig

eig() fungerar endast på reella symmetriska matriser. Icke-symmetriska matriser ger upphov till ValueError. Den returnerar en 2-tupel (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)

Anmärkningar:

  • Egenvärdena kommer tillbaka i ingen särskild ordning. Tillämpa sort() (och samma permutation på egenvektorerna via argsort()) när en sorterad ordning spelar roll.

  • En egenvektor är unik endast upp till en nollskild skalär, så tecknet hos enskilda egenvektorer är inte entydigt definierat. Två korrekta körningar kan producera vektorer med motsatt tecken; detta är ofarligt.

6.13.6. norm

Den euklidiska (Frobenius) normen för en vektor eller matris:

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

Det valfria nyckelordet axis= tar normen längs en enda axel istället för över hela arrayen.

6.13.7. qr

qr() faktoriserar en rektangulär matris A (form (M, N)) till en ortonormal Q och en övre triangulär R så att 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)

Dekompositionen är implementerad via successiva Givens-rotationer. Rätt val för minstakvadratproblem där matrisen inte är symmetrisk.

6.13.8. Lösa system

De två dedikerade lösarna under ulab.scipy.linalg är båda snabbare och mer exakta än np.dot(np.linalg.inv(A), b):

  • solve_triangular(a, b, lower=False)() – lös a @ x = b med antagandet att a är triangulär:

    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)() – givet en Cholesky-faktor L, lös A @ x = b där A = L @ L.T

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

Ta till dessa istället för att invertera närhelst strukturen hos A tillåter – de sparar både elimineringsarbete och den explicita inversen.

6.13.9. Lösa ett litet linjärt system

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

Samma problem kan uttryckas genom att faktorisera A och anropa den lämpliga lösaren – snabbare och mer exakt när A har rätt struktur.

För den fullständiga referensen på argumentnivå, se numpy.linalg — Rutiner för linjär algebra och scipy.linalg — Linjäralgebrarutiner.