6.13. Lineáris algebra

A lineáris algebra egy kamerán kis mátrixokkal végzett munka: egy 3x3-as forgatás, amely egy IMU-mintát a világkoordináta-rendszerbe illeszt, egy kalibrációs mátrix, amely egy lencsét helyrehoz, egy Kalman-szűrő állapot-kovariancia frissítése, egy polinomiális illesztés, amelynek normálegyenletei apró lineáris megoldásként jönnek ki. A numpy.linalg és scipy.linalg almodulok pontosan ezt a méretet fedik le.

A függvények két modulban élnek. A mátrixszorzás-műveletek a numpy legfelső szintjén találhatók; a felbontások és a mátrixinverz a numpy.linalg alatt; a dedikált lineárisrendszer-megoldók a scipy.linalg alatt. Egy 2x2-es mátrixszorzás például:

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. Mi áll rendelkezésre

  • dot() – mátrix- vagy vektorszorzat.

  • cross() – 3D vektor keresztszorzata.

  • trace() – az átló összege.

  • inv() – mátrixinverz.

  • det() – determináns.

  • cholesky() – Cholesky-felbontás (szimmetrikus pozitív definit bemenet).

  • eig() – egy valós szimmetrikus mátrix sajátértékei és sajátvektorai.

  • norm() – egy vektor vagy mátrix 2-normája.

  • qr() – QR-felbontás mode='reduced' (alapértelmezett) vagy mode='complete' móddal.

  • solve_triangular() – az A @ x = b megoldása, ha A háromszögmátrix.

  • cho_solve() – az A @ x = b megoldása az A Cholesky-faktora alapján.

6.13.2. dot, cross, trace

A dot() az, ahogyan a mátrixszorzást a kamerán le kell írni. A @ operátor, amelyet az asztali numpy ugyanerre a feladatra biztosít, nincs megvalósítva a ndarray osztályon, így minden mátrix-vektor, mátrix-mátrix és vektor-skalárszorzat a dot() híváson keresztül megy:

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

A dot() eredménye mindig float dtype-ú. A cross() két 3-vektor keresztszorzata, a trace() pedig egy négyzetes mátrix főátlójának összege.

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

Az inverzet Gauss-Jordan-eliminációval számítja ki, így a inv() ValueError kivételt vált ki, ha a mátrix szinguláris (egy átlóelem nullává válik az elimináció során). A RAM-költség nagyjából a bemenet méretének kétszerese.

A determináns ugyanazt az eliminációt használja újra – a futási idő lényegében megegyezik az inverzével.

Amikor az alkalmazás célja egy lineáris rendszer megoldása, ne invertálj és szorozz – inkább az alábbi dedikált megoldókat részesítsd előnyben. Mindkettő gyorsabb és numerikusan jobban viselkedik.

6.13.4. cholesky

Egy A szimmetrikus pozitív definit mátrixra a cholesky() egy alsó háromszögmátrix L értéket ad vissza úgy, hogy A = L @ L.T

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

Ha a bemenet nem pozitív definit vagy nem szimmetrikus, ValueError kivétel keletkezik.

A Cholesky-faktor egy LU-faktorizáció munkájának fele, és a megfelelő kiindulópont minden olyan problémához, amelynek mátrixáról ismert, hogy szimmetrikus pozitív definit (kovariancia-frissítések, normálegyenletek egy legkisebb négyzetes illesztésből).

6.13.5. eig

A eig() csak valós szimmetrikus mátrixokon működik. A nem szimmetrikus mátrixok ValueError kivételt váltanak ki. Egy 2 elemű (eigenvalues, eigenvectors) rendezett párt ad vissza:

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)

Megjegyzések:

  • A sajátértékek nem meghatározott sorrendben jönnek vissza. Alkalmazd a sort() függvényt (és ugyanazt a permutációt a sajátvektorokra a argsort() segítségével), amikor a rendezett sorrend számít.

  • Egy sajátvektor csak egy nem nulla skalár erejéig egyértelmű, így az egyes sajátvektorok előjele nincs egyértelműen definiálva. Két helyes futtatás ellentétes előjelű vektorokat is produkálhat; ez ártalmatlan.

6.13.6. norm

Egy vektor vagy mátrix euklideszi (Frobenius-) normája:

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

Az opcionális axis= kulcsszó a normát egyetlen tengely mentén veszi, a teljes tömb helyett.

6.13.7. qr

A qr() egy (M, N) alakú A téglalap alakú mátrixot egy ortonormált Q és egy felső háromszögmátrix R szorzatára bont fel úgy, hogy 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)

A felbontás egymást követő Givens-forgatásokon keresztül van megvalósítva. A megfelelő választás olyan legkisebb négyzetes problémákhoz, ahol a mátrix nem szimmetrikus.

6.13.8. Rendszerek megoldása

A két dedikált megoldó az ulab.scipy.linalg alatt egyaránt gyorsabb és pontosabb, mint az np.dot(np.linalg.inv(A), b):

  • solve_triangular(a, b, lower=False)() – az a @ x = b megoldása, feltételezve, hogy a háromszögmátrix:

    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)() – egy L Cholesky-faktor alapján az A @ x = b megoldása, ahol A = L @ L.T

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

Ezeket vedd elő invertálás helyett, amikor A szerkezete megengedi – megtakarítják az eliminációs munkát és az explicit inverzet.

6.13.9. Egy kis lineáris rendszer megoldása

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

Ugyanaz a probléma kifejezhető A faktorizálásával és a megfelelő megoldó hívásával – gyorsabb és pontosabb, ha A rendelkezik a megfelelő szerkezettel.

A teljes argumentumszintű referenciáért lásd: numpy.linalg — Lineáris algebrai rutinok és scipy.linalg — Lineáris algebrai rutinok.