6.13. Linearna algebra#

Linearna algebra na kameri rad je s malim matricama: 3x3 rotacija koja stapa IMU uzorak u svjetski okvir, kalibracijska matrica koja ispravlja leću, ažuriranje kovarijance stanja Kalmanova filtra, polinomijalna prilagodba čije normalne jednadžbe izlaze kao maleno linearno rješavanje. Podmoduli numpy.linalg i scipy.linalg pokrivaju upravo tu razinu.

Funkcije se nalaze u dvama modulima. Operacije matričnog umnoška na su najvišoj razini modula numpy; dekompozicije i inverz matrice nalaze se pod numpy.linalg; namjenski rješavači linearnih sustava nalaze se pod scipy.linalg. Množenje 2-sa-2 matrice, primjerice:

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. Što je dostupno#

  • dot() – matrični ili vektorski umnožak.

  • cross() – vektorski (kros) umnožak u 3-D.

  • trace() – suma dijagonale.

  • inv() – inverz matrice.

  • det() – determinanta.

  • cholesky() – Choleskyjeva dekompozicija (simetričan pozitivno definitan ulaz).

  • eig() – svojstvene vrijednosti i svojstveni vektori realne simetrične matrice.

  • norm() – 2-norma vektora ili matrice.

  • qr() – QR dekompozicija s mode='reduced' (zadano) ili mode='complete'.

  • solve_triangular() – rješava A @ x = b kada je A trokutasta.

  • cho_solve() – rješava A @ x = b uz zadani Choleskyjev faktor matrice A.

6.13.2. dot, cross, trace#

dot() način je na koji se na kameri zapisuje matrično množenje. Operator @ koji stolni numpy pruža za isti posao nije implementiran na ndarray, pa svaki matrično-vektorski, matrično-matrični i vektorski skalarni umnožak ide kroz poziv 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

Rezultat poziva dot() uvijek je dtype-a float. cross() vektorski je umnožak dvaju 3-vektora, a trace() suma glavne dijagonale kvadratne matrice.

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

Inverz se računa Gauss-Jordanovom eliminacijom, pa inv() podiže ValueError kada je matrica singularna (dijagonalni unos postane nula tijekom eliminacije). Trošak RAM-a otprilike je dvostruka veličina ulaza.

Determinanta ponovno koristi istu eliminaciju – vrijeme izvođenja u biti je isto kao za inverz.

Kada je cilj aplikacije riješiti linearni sustav, nemojte invertirati i množiti – radije upotrijebite namjenske rješavače u nastavku. Oba su brža i numerički se bolje ponašaju.

6.13.4. cholesky#

Za simetričnu pozitivno definitnu matricu A, cholesky() vraća donje-trokutastu L takvu da je A = L @ L.T

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

Ako ulaz nije pozitivno definitan ili nije simetričan, podiže se ValueError.

Choleskyjev faktor upola je posla LU faktorizacije i prava je polazna točka za svaki problem čija je matrica poznato simetrična pozitivno definitna (ažuriranja kovarijance, normalne jednadžbe iz prilagodbe metodom najmanjih kvadrata).

6.13.5. eig#

eig() radi samo na realnim simetričnim matricama. Nesimetrične matrice podižu ValueError. Vraća 2-torku (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)

Napomene:

  • Svojstvene vrijednosti vraćaju se bez određenog redoslijeda. Primijenite sort() (te istu permutaciju na svojstvene vektore putem argsort()) kada je sortiran redoslijed bitan.

  • Svojstveni vektor jedinstven je samo do nenul skalara, pa predznak pojedinačnih svojstvenih vektora nije jedinstveno određen. Dva ispravna izvođenja mogu proizvesti vektore suprotnog predznaka; to je bezopasno.

6.13.6. norm#

Euklidska (Frobeniusova) norma vektora ili 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...

Neobavezna ključna riječ axis= uzima normu duž jedne osi umjesto nad cijelim poljem.

6.13.7. qr#

qr() faktorizira pravokutnu matricu A (oblik (M, N)) u ortonormalnu Q i gornje-trokutastu R takve da je 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)

Dekompozicija je implementirana putem uzastopnih Givensovih rotacija. Pravi izbor za probleme najmanjih kvadrata gdje matrica nije simetrična.

6.13.8. Rješavanje sustava#

Dva namjenska rješavača pod ulab.scipy.linalg i brža su i točnija od np.dot(np.linalg.inv(A), b):

  • solve_triangular(a, b, lower=False)() – rješava a @ x = b pretpostavljajući da je a trokutasta:

    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)() – uz zadani Choleskyjev faktor L, rješava A @ x = b gdje je A = L @ L.T

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

Posegnite za ovima umjesto invertiranja kad god struktura matrice A to dopušta – štede posao eliminacije i eksplicitni inverz.

6.13.9. Rješavanje malog linearnog sustava#

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

Isti se problem može izraziti faktorizacijom matrice A i pozivom odgovarajućeg rješavača – brže i točnije kada A ima pravu strukturu.

Za potpunu referencu na razini argumenata pogledajte numpy.linalg — Rutine linearne algebre i scipy.linalg — Rutine linearne algebre.