Rozwiązywanie układów równań liniowych

Układ równań liniowych można zapisać w postaci macierzowej jako Ax = b, gdzie A jest macierzą współczynników, x jest wektorem niewiadomych, a b jest wektorem prawej strony. Na przykład dwa równania z dwiema niewiadomymi można wyrazić jako równanie macierzowe z macierzą 2×2. Rozwiązanie układu oznacza znalezienie x takiego, że A pomnożone przez x jest równe b. Funkcja np.linalg.solve(A, b) w NumPy wykonuje to wydajnie i dokładnie.

Matematycznie rozwiązaniem równania Ax = b jest x = A⁻¹b. Jednak obliczenie np.linalg.inv(A) @ b jest wolniejsze i mniej dokładne niż bezpośrednie wywołanie np.linalg.solve(A, b). Funkcja solve() wewnętrznie używa rozkładu LU, który wymaga mniejszej liczby operacji i powoduje kumulowanie mniejszych błędów zmiennoprzecinkowych. Zasada praktyczna: jeśli rozwiązują Państwo równanie Ax=b względem x, należy zawsze używać solve, a nigdy inv.

np.linalg.solve(A, b) wymaga, aby A była kwadratowa i nieosobliwa. Jeśli A jest osobliwa (wyznacznik wynosi zero), funkcja zgłasza np.linalg.LinAlgError: Singular matrix. Macierz A jest osobliwa, gdy jej wiersze są liniowo zależne — na przykład gdy jedno równanie jest po prostu wielokrotnością innego. Przed wywołaniem solve należy zawsze sprawdzić, czy układ jest poprawnie określony. Jeśli A nie jest kwadratowa, należy zamiast tego użyć np.linalg.lstsq().

np.linalg.solve wewnętrznie stosuje rozkład LU: rozkłada A na macierz trójkątną dolną L i macierz trójkątną górną U, tak że A = LU. Rozwiązanie sprowadza się wtedy do dwóch prostszych układów trójkątnych (podstawiania w przód i wstecz), które są znacznie wydajniejsze niż wielokrotne wykonywanie eliminacji Gaussa od początku. Jest to ten sam algorytm, którego używają MATLAB, Fortran LAPACK i praktycznie każda biblioteka do obliczeń naukowych.

np.linalg.solve(A, B) może obsługiwać macierz B z wieloma kolumnami, rozwiązując układ jednocześnie dla każdej kolumny. Jeśli B ma kształt (n, k), wynik X ma kształt (n, k), przy czym każda kolumna X jest rozwiązaniem odpowiadającym właściwej kolumnie B. Jest to przydatne w symulacjach fizycznych i uczeniu maszynowym, gdy ten sam układ współczynników trzeba rozwiązać dla wielu różnych prawych stron.

Gdy układ jest nadokreślony (ma więcej równań niż niewiadomych), nie istnieje rozwiązanie dokładne. np.linalg.lstsq(A, b) znajduje takie x, które minimalizuje sumę kwadratów reszt ||Ax - b||² — dokładnie to oblicza regresja liniowa. Funkcja zwraca rozwiązanie x, reszty, rząd oraz wartości osobliwe. Argument rcond=None wycisza ostrzeżenie o przestarzałym użyciu i ustawia próg, poniżej którego wartości osobliwe są traktowane jako zero.

Równania normalne zwykłej metody najmniejszych kwadratów mają postać A.T @ A @ x = A.T @ b, czyli kwadratowego układu, który można rozwiązać za pomocą np.linalg.solve. Daje to taki sam wynik jak lstsq, gdy A.T @ A jest dobrze uwarunkowana. W praktyce preferowane jest lstsq, ponieważ używa SVD i pozostaje stabilne nawet wtedy, gdy kolumny są prawie współliniowe, ale równania normalne pomagają zrozumieć matematykę stojącą za regresją.

Przed rozwiązaniem warto sprawdzić współczynnik uwarunkowania za pomocą np.linalg.cond(A). Współczynnik bliski 1 oznacza, że układ jest dobrze uwarunkowany, a rozwiązanie — wiarygodne. Bardzo duży współczynnik (np. 1e12) oznacza, że macierz jest prawie osobliwa: małe zmiany b powodują ogromne zmiany x, a rozwiązanie może być numerycznie niewiarygodne. W takim przypadku należy rozważyć regularyzację (regresję ridge) albo użycie bardziej odpornego solvera.

W przypadku bardzo dużych układów, w których A ma głównie elementy zerowe (jest macierzą rzadką), przechowywanie pełnej macierzy marnuje pamięć, a użycie gęstej funkcji np.linalg.solve jest powolne. SciPy udostępnia scipy.sparse oraz scipy.sparse.linalg.spsolve, które wykorzystują rzadkość macierzy w ogromnych układach, takich jak modele elementów skończonych lub problemy sieciowe. Rozpoznanie, kiedy macierz jest rzadka, jest ważnym aspektem wydajności w obliczeniach naukowych.

Po uzyskaniu wektora rozwiązania x za pomocą np.linalg.solve(A, b) należy zawsze zweryfikować rozwiązanie, obliczając A @ x i porównując wynik z b za pomocą np.allclose(A @ x, b). W arytmetyce zmiennoprzecinkowej dokładna równość występuje rzadko — allclose sprawdza zgodność w ramach niewielkiej tolerancji. Norma residuum ||Ax - b|| powinna być bliska iloczynowi epsilon maszynowego i skali problemu dla dobrze uwarunkowanego układu.

Solvery liniowe są wykorzystywane w całym uczeniu maszynowym: regresja liniowa rozwiązuje równania normalne, procesy Gaussa rozwiązują układy z macierzami jąder, a filtry Kalmana rozwiązują równania aktualizacji stanu. Rozumiejąc działanie np.linalg.solve i lstsq, mogą Państwo implementować te algorytmy od podstaw i diagnozować problemy, gdy wyniki są nieoczekiwane. Moduł linalg biblioteki NumPy udostępnia te same narzędzia numeryczne, które są używane w produkcyjnych bibliotekach naukowych.

Sprawdź swoją wiedzę na temat koncepcji analizy danych z tej lekcji.

W tej lekcji nauczyli się Państwo, że: np.linalg.solve(A, b) efektywnie rozwiązuje kwadratowe układy równań liniowych za pomocą rozkładu LU, np.linalg.lstsq() znajduje rozwiązanie metodą najmniejszych kwadratów dla układów nadokreślonych, a liczba uwarunkowania mierzy stabilność numeryczną — duże wartości ostrzegają o macierzach bliskich osobliwym. W następnej części omówimy wartości własne i SVD, które stanowią podstawę PCA oraz wielu algorytmów uczenia maszynowego.