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.