6.13. Algebra liniowa

Algebra liniowa na kamerze to praca z małymi macierzami: obrót 3x3 łączący próbkę IMU z układem świata, macierz kalibracji prostująca obiektyw, aktualizacja kowariancji stanu w filtrze Kalmana, dopasowanie wielomianu, którego równania normalne sprowadzają się do drobnego rozwiązania liniowego. Podmoduły numpy.linalg i scipy.linalg obejmują dokładnie tę skalę.

Funkcje znajdują się w dwóch modułach. Operacje iloczynu macierzowego są na najwyższym poziomie numpy; rozkłady i odwrotność macierzy znajdują się pod numpy.linalg; dedykowane solvery układów liniowych znajdują się pod scipy.linalg. Na przykład mnożenie macierzy 2 na 2:

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. Co jest dostępne

  • dot() – iloczyn macierzy lub wektorów.

  • cross() – iloczyn wektorowy w 3-W.

  • trace() – suma przekątnej.

  • inv() – odwrotność macierzy.

  • det() – wyznacznik.

  • cholesky() – rozkład Choleskiego (wejście symetryczne dodatnio określone).

  • eig() – wartości i wektory własne rzeczywistej macierzy symetrycznej.

  • norm() – norma 2 wektora lub macierzy.

  • qr() – rozkład QR z mode='reduced' (domyślnie) lub mode='complete'.

  • solve_triangular() – rozwiązuje A @ x = b, gdy A jest trójkątna.

  • cho_solve() – rozwiązuje A @ x = b dla danego czynnika Choleskiego A.

6.13.2. dot, cross, trace

dot() to sposób zapisu mnożenia macierzy na kamerze. Operator @, który komputerowy numpy udostępnia do tego samego zadania, nie jest zaimplementowany na ndarray, więc każdy iloczyn macierz-wektor, macierz-macierz oraz iloczyn skalarny wektorów przechodzi przez wywołanie 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

Wynik dot() jest zawsze typu dtype float. cross() to iloczyn wektorowy dwóch 3-wektorów, a trace() to suma głównej przekątnej macierzy kwadratowej.

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

Odwrotność jest obliczana metodą eliminacji Gaussa-Jordana, więc inv() zgłasza ValueError, gdy macierz jest osobliwa (podczas eliminacji wpis na przekątnej staje się zerem). Koszt pamięci RAM to mniej więcej dwukrotność rozmiaru wejścia.

Wyznacznik ponownie wykorzystuje tę samą eliminację – czas działania jest zasadniczo taki sam jak dla odwrotności.

Gdy celem aplikacji jest rozwiązanie układu liniowego, nie odwracaj i nie mnóż – preferuj dedykowane solvery poniżej. Oba są szybsze i lepiej zachowują się numerycznie.

6.13.4. cholesky

Dla symetrycznej dodatnio określonej macierzy A funkcja cholesky() zwraca dolnotrójkątną L taką, że A = L @ L.T

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

Jeśli wejście nie jest dodatnio określone lub nie jest symetryczne, zgłaszany jest ValueError.

Czynnik Choleskiego to połowa pracy rozkładu LU i jest właściwym punktem wyjścia dla każdego problemu, którego macierz jest znana jako symetryczna dodatnio określona (aktualizacje kowariancji, równania normalne z dopasowania metodą najmniejszych kwadratów).

6.13.5. eig

eig() działa tylko na macierzach rzeczywistych symetrycznych. Macierze niesymetryczne zgłaszają ValueError. Zwraca 2-elementową krotkę (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)

Uwagi:

  • Wartości własne zwracane są w nieokreślonej kolejności. Zastosuj sort() (oraz tę samą permutację do wektorów własnych za pomocą argsort()), gdy posortowana kolejność ma znaczenie.

  • Wektor własny jest jednoznaczny tylko z dokładnością do niezerowego skalara, więc znak poszczególnych wektorów własnych nie jest jednoznacznie określony. Dwa poprawne uruchomienia mogą dać wektory o przeciwnym znaku; jest to nieszkodliwe.

6.13.6. norm

Norma euklidesowa (Frobeniusa) wektora lub macierzy:

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

Opcjonalne słowo kluczowe axis= oblicza normę wzdłuż pojedynczej osi zamiast po całej tablicy.

6.13.7. qr

qr() rozkłada prostokątną macierz A (kształt (M, N)) na ortonormalną Q i górnotrójkątną R takie, że 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)

Rozkład jest zaimplementowany za pomocą kolejnych obrotów Givensa. Właściwy wybór dla problemów najmniejszych kwadratów, gdy macierz nie jest symetryczna.

6.13.8. Rozwiązywanie układów

Dwa dedykowane solvery pod ulab.scipy.linalg są zarówno szybsze, jak i dokładniejsze niż np.dot(np.linalg.inv(A), b):

  • solve_triangular(a, b, lower=False)() – rozwiązuje a @ x = b przy założeniu, że a jest trójkątna:

    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)() – dla danego czynnika Choleskiego L rozwiązuje A @ x = b, gdzie A = L @ L.T

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

Sięgaj po nie zamiast odwracać, gdy tylko pozwala na to struktura A – oszczędzają pracę eliminacji oraz jawnej odwrotności.

6.13.9. Rozwiązywanie małego układu liniowego

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

Ten sam problem można wyrazić, rozkładając A i wywołując odpowiedni solver – szybciej i dokładniej, gdy A ma właściwą strukturę.

Pełną dokumentację na poziomie argumentów znajdziesz w numpy.linalg — Procedury algebry liniowej oraz scipy.linalg — Procedury algebry liniowej.