Złożoność inwersji macierzy w liczbach

Rozwiązuję równania różniczkowe, które wymagają odwrócenia gęstych macierzy kwadratowych. Ta inwersja macierzy zużywa najwięcej mojego czasu obliczeniowego, więc zastanawiałem się, czy używam najszybszego dostępnego algorytmu.

Mój obecny wybór to numpy.linalg.inv . Z moich danych liczbowych wynika, że skaluje się jako gdzie n jest liczbą rzędów, więc metoda wydaje się być eliminacją Gaussa. $O(n^3)$

Według Wikipedii dostępne są szybsze algorytmy. Czy ktoś wie, czy istnieje biblioteka, która je implementuje?

Zastanawiam się, dlaczego numpy nie używa tych szybszych algorytmów?

numpy dense-matrix inverse

— fizyka
źródło

Musisz wcześniej wykonać swoje macierze. Spójrz na Scipy. Rzadko za twoją pomoc. Zawiera wiele potrzebnych narzędzi.

— Tobal

@Tobal nie jestem pewien, czy śledzę ... jak „wykonałbyś” matrycę? i dokładnie jak by to scipy.sparsepomogło?

— GoHokies

@GoHokies scipy jest uzupełnieniem numpy. Matryce gęste / rzadkie muszą zostać zaimplementowane na długo przed wykonaniem niektórych obliczeń, co usprawnia obliczenia. Proszę przeczytać ten docs.scipy.org/doc/scipy/reference/sparse.html, który najlepiej wyjaśnia niż mój zły angielski.

— Tobal

@Tobal Pytanie odnosi się konkretnie do gęstych matryc, więc nie rozumiem, jak scipy.sparseto jest istotne?

— Christian Clason

@Tobal - Myślę, że nadal nie rozumiem. Co dokładnie masz na myśli przez „preformuj swoje macierze” i „macierze muszą zostać zaimplementowane na długo przed wykonaniem obliczeń”? Jeśli chodzi o twój ostatni komentarz, na pewno zgodzisz się, że techniki, które można zastosować do rzadkich i gęstych matryc, są bardzo różne.

— Wolfgang Bangerth,

Odpowiedzi:

(Robi się zbyt długo na komentarze ...)

Zakładam, że faktycznie musisz obliczyć odwrotność w swoim algorytmie. ¹ Po pierwsze, należy zauważyć, że te alternatywne algorytmy nie są tak naprawdę twierdzone, że są szybsze , tylko że mają lepszą asymptotyczną złożoność (co oznacza, że wymagana liczba elementarnych operacji rośnie wolniej). W rzeczywistości są one w rzeczywistości (znacznie) wolniejsze niż standardowe podejście (dla danego ) z następujących powodów: $n$

-notation skóry stałym przed mocy , które mogą być astronomically duża - tak duża, że może być znacznie mniejszy niż dla każdego , że może być obsługiwany przez dowolny komputer w dającej się przewidzieć przyszłości. (Tak jest na przykład w przypadku algorytmu Coppersmith – Winograd.) $\mathcal{O}$ $n$ $C_1 n^3$ $C_2 n^{2.x}$ $n$
Złożoność zakłada, że każda (arytmetyczna) operacja zajmuje ten sam czas - ale w praktyce nie jest to prawdą: Mnożenie wiązki liczb o tej samej liczbie jest znacznie szybsze niż pomnożenie tej samej liczby różnych liczb. Wynika to z faktu, że głównym wąskim gardłem w bieżącym przetwarzaniu danych jest umieszczanie danych w pamięci podręcznej, a nie faktyczne operacje arytmetyczne na tych danych. Dlatego algorytm, który można zmienić tak, aby miał pierwszą sytuację (zwaną pamięcią podręczną ), będzie znacznie szybszy niż ten, w którym nie jest to możliwe. (Tak jest na przykład w przypadku algorytmu Strassena).

Ponadto stabilność numeryczna jest co najmniej tak samo ważna jak wydajność; i tutaj znowu standardowe podejście zwykle wygrywa.

Z tego powodu standardowe wysokowydajne biblioteki (BLAS / LAPACK, które Numpy wywołuje, gdy poprosi go o obliczenie odwrotności) zwykle implementują to podejście. Oczywiście istnieją implementacje Numpy np. Algorytmu Strassena, ale algorytm ręcznie dostrojony na poziomie asemblera wyraźnie pokona algorytm napisany w języku wysokiego poziomu dla dowolnego rozsądnego rozmiaru matrycy. $\mathcal{O}(n^3)$ $\mathcal{O}(n^{2.x})$

¹ Byłbym jednak zły, gdybym nie zauważył, że jest to bardzo rzadko naprawdę konieczne: za każdym razem, gdy trzeba obliczyć produkt , należy zamiast tego rozwiązać układ liniowy (np. using ) i zamiast tego użyj - jest to znacznie bardziej stabilne i można to zrobić (w zależności od struktury macierzy ) znacznie szybciej. Jeśli musisz użyć wiele razy, możesz wstępnie obliczyć faktoryzację (która jest zwykle najdroższą częścią rozwiązania) i użyć go później.

A^{- 1} b

$A^{-1}b$

A x = b

$Ax=b$ numpy.linalg.solve

x

$x$

A

$A$

A^{- 1}

$A^{-1}$

A

$A$

— Christian Clason
źródło

Świetna odpowiedź, dziękuję panu, w szczególności za wskazanie diabła w szczegółach (stałe w dużej notacji O), która robi dużą różnicę między prędkością teoretyczną a prędkością praktyczną.

— gaboryczny

Myślę, że część „odwrotność jest rzadko konieczna” powinna zostać podkreślona bardziej. Jeśli celem jest rozwiązanie układu równań różniczkowych, nie wydaje się prawdopodobne, że potrzebna jest pełna odwrotność.

— Jared Goguen,

@ o_o Cóż, to był mój pierwszy oryginalny komentarz (który usunąłem po scaleniu ich wszystkich w jedną odpowiedź). Ale pomyślałem, że z korzyścią dla strony (i późniejszych czytelników) odpowiedź powinna odpowiedzieć na rzeczywiste pytanie w pytaniu (które jest zarówno uzasadnione, jak i tematyczne), nawet jeśli kryje się za tym problem XY. Poza tym nie chciałem zabrzmieć zbyt upomniająco ...

— Christian Clason,

n

$n$

A

$A$

Prawdopodobnie powinieneś zauważyć, że pochowany głęboko w kodzie źródłowym numpy (patrz https://github.com/numpy/numpy/blob/master/numpy/linalg/umath_linalg.c.src ) procedura inv próbuje wywołać funkcję dgetrf z systemowego pakietu LAPACK, który następnie dokonuje dekompozycji LU oryginalnej macierzy. Jest to moralnie równoważne z eliminacją Gaussa, ale można je dostroić do nieco mniejszej złożoności, stosując szybsze algorytmy mnożenia macierzy w wysokowydajnym BLAS.

Jeśli podążysz tą drogą, powinieneś zostać ostrzeżony, że zmuszanie całego łańcucha bibliotek do korzystania z nowej biblioteki zamiast systemowej dostarczonej z twoją dystrybucją jest dość złożone. Jedną alternatywą we współczesnych systemach komputerowych jest spojrzenie na równoległe metody przy użyciu pakietów takich jak scaLAPACK lub (w świecie python) petsc4py. Jednak zazwyczaj są one szczęśliwsze, gdy są stosowane jako solwery iteracyjne dla systemów algebry liniowej, niż stosowane do metod bezpośrednich, a PETSc w szczególności do układów rzadkich bardziej niż gęstych.

— origimbo
źródło