Jak zmienić kolejność zmiennych, aby uzyskać pasmową macierz minimalnej przepustowości?

15

Próbuję rozwiązać równanie Poissona 2D na podstawie różnic skończonych. W tym procesie otrzymuję rzadką macierz z tylko zmiennymi w każdym równaniu. Na przykład, jeśli zmienne byłyby , dyskretyzacja dałaby: $5$ $U$

U_{i - 1, j} + U_{i + 1, j} - 4 U_{i, j} + U_{i, j - 1} + U_{i, j + 1} = f_{i, j}

$U_{i-1,j} + U_{i+1,j} -4U_{i,j} + U_{i,j-1} + U_{i,j+1} = f_{i,j}$

Wiem, że mogę rozwiązać ten system za pomocą metody iteracyjnej, ale przyszło mi do głowy, że jeśli odpowiednio uporządkuję zmienne, być może uda mi się uzyskać macierz pasmową, którą można rozwiązać metodą bezpośrednią (tj. Eliminacja Gaussa w / o pivoting). czy to możliwe? Czy istnieją jakieś strategie pozwalające to zrobić w przypadku innych, być może mniej ustrukturyzowanych rzadkich systemów?

— Paweł
źródło

2

Coś takiego jak Cuthill-McKee?

— JM

Ciekawe ... Nigdy wcześniej nie słyszałem o algorytmie Cuthill-McKee! :)

— Paweł

1

Istnieje również Reverse Cuthill-McKee.

— Geoff Oxberry

1

Mam nadzieję, że wynika to z odpowiedzi, ale nie chcesz używać rozwiązania pasmowego dla tego problemu ani wybierać kolejności, która minimalizuje przepustowość. Być może pytanie lub wybraną odpowiedź można zredagować, aby to wyjaśnić, w przeciwnym razie obawiam się, że ten mit zostanie utrwalony. Dałem porównanie wizualne i porównałem wypełnij scicomp.stackexchange.com/a/880/119 .

— Jed Brown

@JedBrown: Właściwie nie całkiem pracuję z problemem Poissona per se ... Mój problem ma podobną strukturę do problemu Poissona ... Wskaźniki zmiennych (i's i j's) są dokładnie takie same i matryca jest dominująca po przekątnej, a wpisy poza przekątną (w tym samym rzędzie) dodają dokładnie do sumy wejścia po przekątnej.

— Paweł

13

Jest to dobrze zbadany problem w dziedzinie solverów rzadkich. Gorąco polecam zapoznanie się z omówieniem metody wielopłaszczyznowej przez Josepha Liu , aby lepiej zrozumieć, w jaki sposób zmiany kolejności i supernody wpływają na czas wypełniania i rozwiązania.

Zagnieżdżone rozwarstwienie jest niezwykle powszechnym sposobem generowania zmiany kolejności i zasadniczo polega na rekurencyjnym dzieleniu grafów. MeTiS jest de facto standardem do partycjonowania grafów i możesz przeczytać o niektórych pomysłach tutaj . Innym często używanym pakietem jest SCOTCH , a Chaco jest również ważne, ponieważ jego autorzy wprowadzili wielopoziomowe partycjonowanie grafów , co jest również podstawową ideą MeTiS .

George i Liu wykazały ich klasycznym książki że 2d rozwiązania rzadki direct wymagają jedynie pracy i pamięci, natomiast 3d rzadki direct wymaga praca i pamięć . $O(n^{3/2})$ $O(n \log n)$ $O(n^2)$ $O(n^{4/3})$

— Jack Poulson
źródło

Czy masz cytat z referencji George'a i Liu?

— Paweł

Dodany; Już miałam wysiąść z samochodu, kiedy go po raz pierwszy przesłałam. Wiem, że gdzieś istnieje darmowa wersja książki (Jed wie, gdzie ona jest), ale nie mogłem jej znaleźć.

— Jack Poulson,

Zaktualizowałem link, aby wskazywał plik PDF książki zamiast recenzji książki.

— Jed Brown

@JedBrown To była świetna referencja! Dzięki wielkie! :)

— Paul

1

@Alexander Wszyscy przypisują 3D związaną z George'em i Liu, choć nie wiem, czy wyraźnie to zaznaczają w książce. Jest to jednak oczywiste z teorii. Minimalna separatora wierzchołek dla

siatki jest

. Gęsta matryca wiąże się z tym Supernode ma

dane i wymaga

n = m \times m \times m

$n = m\times m\times m$

n^{2 / 3} = m \times m

$n^{2/3} = m\times m$

(n^{2 / 3})^{2} = n^{4 / 3}

$(n^{2/3})^2 = n^{4/3}$

(n^{2 / 3})^{3} = n^{2}

$(n^{2/3})^3 = n^2$ operacje do uwzględnienia. Pojęcie logarytmiczne w przypadku 2D jest bardziej subtelne i jest omówione w rozdziale 8 dotyczącym zagnieżdżenia, w którym osiągnięto dolną granicę.

— Jed Brown

5

Cuthill-McKee to de facto standard tego, co chcesz robić. Jeśli chcesz bawić się tą metodą, istnieje łatwa w użyciu implementacja algorytmu (i jego odwrotności) w bibliotece grafów Boost (BGL), a dokumentacja zawiera przykłady, jak go używać.

— Wolfgang Bangerth
źródło

faktycznie odwraca Cuhill-McKee; zwykle daje mniej wypełnienia. Ale zagnieżdżona kolejność sekcji jest znacznie lepsza niż kolejność niskiej przepustowości.

— Arnold Neumaier

4

Mówiąc o metodach wielopłaszczyznowych, Tim Davis , który pracuje nad wielopłaszczyznowymi metodami faktoryzacji LU ( UMFPACK ), ma szereg procedur, które będą zmieniać kolejność matryc w celu zminimalizowania wypełnienia. Można je znaleźć tutaj, jako część SuiteSparse . SuiteSparse korzysta z MeTiS.

Jeszcze jedna rzecz do zapamiętania: w niektórych problemach możesz być sprytny w porządkowaniu zmiennych, aby uzyskać pasmowe lub zbliżone do pasmowych wzorce, które mogą zaoszczędzić ci kłopotów (i czasu procesora) wywołania tych algorytmów. Jednak ta sprytna zmiana kolejności wymaga wglądu z twojej strony i nie jest tak ogólna, jak algorytmy zmiany kolejności oparte na teorii grafów, o których wspominali ludzie w swoich odpowiedziach tutaj.

— Geoff Oxberry
źródło

Nie ma za co, Paul. Jeśli ci się podoba, zagłosuj.

— Geoff Oxberry

3

W zastosowanych kręgach matematycznych jest algorytm o nazwie ADI (Alternative Direction Implicit) i operator Split w kręgach fizyki, który robi właściwie to, co opisujesz. Jest to metoda iteracyjna i postępuje zgodnie z tą podstawową procedurą:

$y$ $x$
$x$ $y$
Powtarzaj 1 i 2, aż błąd będzie tak mały, jak chcesz.

Nie znam formalnej złożoności tego algorytmu, ale zauważyłem, że za każdym razem, gdy go używam, zbiega się w mniejszej liczbie iteracji niż rzeczy takie jak Jacobi i Gauss-Seidel.

— Dan
źródło

2

Jeśli zdecydujesz się wybrać trasę podziału operatora, należy zachować ostrożność, ponieważ wiadomo, że podział operatora powoduje błędy w ustalonych rozwiązaniach w niektórych przypadkach. (Jeden z moich kolegów z laboratorium opracował sposób na przezwyciężenie tej trudności, ale nie sądzę, aby ją opublikował.) Wiadomo również, że podział operatora powoduje błędy numeryczne. Istnieją ugruntowane sposoby oszacowania tych błędów a posteriori ; Don Estep wykonał doskonałą pracę w tej dziedzinie.

— Geoff Oxberry

@GeoffOxberry Wygląda na to, że mówisz o innym podziale. Możesz użyć ADI w całkowicie niejawnym schemacie, który nie zawiera błędu podziału, ponieważ faktycznie rozwiązuje system. Istnieją również metody IMEX, które rygorystycznie kontrolują błędy podziału.

— Jed Brown

x

$x$

y

$y$

Nigdy nie słyszałem o dzieleniu Godunowa i Stranga. Zwykle dzielę operatora na wzór Baker-Campbell-Hausdorf. Czy to to samo?

— Dan