Czy rozwiązanie liniowego układu równań można aproksymować tylko dla kilku pierwszych zmiennych?

15

Mam liniowy układ równań wielkości mxm, gdzie m jest duży. Jednak zmienne, które mnie interesują, to tylko pierwsze n zmiennych (n jest małe w porównaniu do m). Czy istnieje sposób, aby zbliżyć rozwiązanie dla pierwszych m wartości bez konieczności rozwiązywania całego systemu? Jeśli tak, czy zbliżenie to byłoby szybsze niż rozwiązanie pełnego układu liniowego?

Paweł
źródło
2
Nie, chyba że twoja funkcja forsowania jest również ograniczona do pierwszych n zmiennych. Jeśli tak, możesz utworzyć dopełnienie Schura, choć prawdopodobnie jest ono gęste. Jeśli twój pierwotny operator jest rzadki, może nie być tego wart.
Jack Poulson,
1
Przypuszczam, że można użyć eliminacji gaussowskiej, zaczynając od prawego dolnego rogu matrycy. Byłoby to ~ 2x szybsze niż zwykła eliminacja gaussa, jeśli zależy ci tylko na kilku pierwszych elementach i zatrzymasz się w połowie. Nie wiem, jak by to się porównało z metodami iteracyjnymi.
Dan.
4
@OscarB: Proszę nie. Reguła Cramera jest okrucieństwem w arytmetyki zmiennoprzecinkowej. Nigdy nie słyszałem o tym, aby był używany do poważnych obliczeń i potrzeba sporo uwagi, aby uniknąć złożoności czynnikowej , gdzie wciąż nie konkuruje on z eliminacją Gaussa.
Jack Poulson,
1
@Paul: Większość redukcji zamówień modeli jest wykorzystywana w kontekście dużych systemów ODE lub DAE. Czasami metodologie redukcji są motywowane przez systemy ODE lub DAE, które wynikają z dyskretyzacji PDE. Nie widziałem redukcji modelu stosowanej w równaniach czysto algebraicznych. (Jeśli masz, prześlij mi referencje, ponieważ robię moją pracę magisterską na temat metod redukcji modeli i byłbym bardzo zainteresowany, aby to zobaczyć.) Jeśli chcesz, mógłbym naszkicować, jak może wyglądać redukcja modelu, jeśli traktujemy równania algebraiczne jako zdegenerowany przypadek układu równań różniczkowo-algebraicznych.
Geoff Oxberry
1
@JackPoulson - czy masz coś przeciwko podsumowaniu komentarza jako odpowiedzi? Myślę, że jest to najbardziej poprawne rozwiązanie i nie chcę, aby zgubiło się w komentarzach.
Aron Ahmadia,

Odpowiedzi:

13

Jak zauważyli inni, jest to trudne w przypadku bezpośredniego solvera. To powiedziawszy, nie jest tak trudno zrobić z iteracyjnymi rozwiązaniami. W tym celu należy zauważyć, że większość iteracyjnych solverów w taki czy inny sposób minimalizuje błąd w odniesieniu do niektórych norm. Często norma ta jest albo indukowana przez samą matrycę, ale czasami jest to również norma wektora l2. Ale nie musi tak być: możesz wybrać normę, w której chcesz zminimalizować błąd (lub resztkowy), i możesz na przykład wybrać normę, w której ważisz komponenty, na których Ci zależy, za pomocą 1 i wszystkie inne z 1e-12, tj. na przykład coś takiego (1e-24)N i = 6 x 2 i i odpowiedni iloczyn skalarny. Następnie napisz wszystkie kroki iteracyjnego solvera w odniesieniu do tej normy i iloczynu skalarnego, a otrzymasz iteracyjny solver, który zwraca znacznie większą uwagę na elementy wektorowe, na których ci zależy, niż na inne.||x||2)=ja=15xja2)+ja=6N.xja2)

Pytanie oczywiście brzmi, czy potrzebujesz mniej iteracji niż w przypadku produktu norm / skalar, który waży wszystkie składniki jednakowo. Ale tak powinno być: powiedzmy, że zależy ci tylko na pięciu pierwszych elementach wektorowych. Następnie powinieneś potrzebować maksymalnie pięciu iteracji, aby zmniejszyć błąd o współczynnik 1e12, ponieważ pięć iteracji jest potrzebne dla opisującego je systemu 5x5. To nie jest dowód, ale jestem całkiem pewien, że rzeczywiście powinieneś uciec o znacznie mniejszej liczbie iteracji, jeśli waga w normie (1e-12 powyżej) jest mniejsza niż tolerancja, z którą chcesz iteracyjnie rozwiązać układ liniowy .

Wolfgang Bangerth
źródło
2
Hmm, dobra racja. Byłbym zainteresowany zobaczeniem prawdziwego przykładu, ponieważ trochę martwię się efektami jedynie próby rozwiązania kilku stopni swobody; nawet jeśli reszta może być mała, być może norma błędu jest wciąż dość duża (zrób to, aby skutecznie zignorować większość operatora).
Jack Poulson
Intuicyjnie wydaje się, że działa to tylko wtedy, gdy elementy bardzo małego systemu naprawdę dominują w odpowiedzi w sensie L2 (lub normie, w której rozumiesz błąd, który chcesz mierzyć). W przeciwnym razie uważam, że troska Jacka jest słuszna, ale na pewno byłbym zainteresowany, aby zobaczyć liczbowy dowód na to ...
Aron Ahmadia,
Trzeba się upewnić, że wybierzesz metodę, która minimalizuje błąd , a nie resztę. Myślę, że MinErr może być dobrym punktem wyjścia.
Wolfgang Bangerth
@WolfgangBangerth: Nie znam MINERR: czy to jest główne odniesienie?
Jack Poulson
1
Nawet to nie wystarczy, ponieważ będziesz niedokładny. Nie można uzyskać kilku składników dokładnie przy użyciu tego ważenia.
Matt Knepley
17

Formowanie uzupełnienia Schur

Załóżmy, że permutowałeś i podzieliłeś swoją macierz na formę

A=(A11A12A21A22),

tak, że zawiera stopnie swobody zainteresowania i jest znacznie mniejszy niż A 11 , wówczas można utworzyć dopełnienie SchurA22A11

S22:=A22A21A111A12,

albo poprzez częściową prawostronną LU faktoryzację, albo wyraźną formułę, a następnie można zrozumieć w następującym znaczeniu:S22

S22x=y(A11A12A21A22)(x)=(0y),

gdzie oznacza „nieciekawą” część rozwiązania. Tak więc, pod warunkiem, że prawa strona, która jest tylko niezerowa w stopniach swobody uzupełnienia Schur S 22 , musimy rozwiązać tylko S 22 , aby uzyskać część rozwiązania odpowiadającą tym stopniom swobody.S22S22

Złożoność obliczeniowa w nieustrukturyzowanym gęstym przypadku

Ustawienie do wysokości A i n do wysokości A 22 , standardowe metody obliczeniowej S 22 jest Czynnikiem l 11 U 11 : = 11 (zignorujmy odchylany do tej pory) w przybliżeniu 2 / 3 ( N - n ) 3 prace, a następnie do formyNAnA22S22L11U11:=A112/3(Nn)3

S22:=A22(A21U111)(L111A12)=A22A21A111A12

używając dwóch trójkątów wymagających pracy , a następnie wykonując aktualizację do A 22 w pracy 2 n 2 ( N - n ) .n(Nn)2A222n2(Nn)

Zatem całkowity pracy jest w przybliżeniu . Gdy n jest bardzo mała, N - n N , więc koszty te mogą być widoczne w przybliżeniu 2 / 3 N 3 , który jest koszt pełnego faktoryzacji.2/3(Nn)3+2n(Nn)2+2n2(Nn)nNnN2/3N3

Korzyścią jest to, że jeśli istnieje bardzo duża liczba prawostronnych stron do rozwiązania za pomocą tego samego układu równań, wówczas może potencjalnie zostać ponownie użyty wiele razy, przy czym każde rozwiązanie wymagałoby tylko 2 n 2 pracy (zamiast pracy 2 N 2 ), jeśli S 22 jest uwzględniony.S222n22N2S22

Złożoność obliczeniowa w (typowym) rzadkim przypadku

Jeśli twój rzadki system powstał z jakiegoś rodzaju przybliżenia skończonej różnicy lub przybliżenia elementu skończonego, wówczas solwery rzadkie bezpośrednie prawie na pewno będą w stanie wykorzystać część struktury; Systemy 2d można rozwiązać pracy i O ( N log N ) podczas przechowywania, natomiast systemy 3D może być rozwiązany O ( N 2 ) pracy i O ( N 4 / 3 ) pamięci. Systemy oparte na faktorach można następnie rozwiązać przy takim samym nakładzie pracy, jak wymagania dotyczące pamięci.O(N3/2)O(NlogN)O(N2)O(N4/3)

Przywołanie złożoności obliczeniowej polega na tym, że jeśli i masz układ 2d, a ponieważ uzupełnienie Schur będzie prawdopodobnie gęste, złożoność rozwiązania przy uwzględnieniu faktoryzowanego uzupełnienia Schur będzie wynosićO(n2)=O(N), w którym brakuje tylko współczynnika logarytmicznego w porównaniu do rozwiązania pełnego system! 3D, wymagaO(N)pracy zamiastO(N 4 / 3 ).nNO(n2)=O(N)O(N)O(N4/3)

Dlatego ważne jest, aby pamiętać, że w twoim przypadku, gdzie , znaczne oszczędności przyniesie tylko praca w kilku wymiarach i rozwiązanie wielu problemów po prawej stronie.n=N

Jack Poulson
źródło
1
To świetne podsumowanie metody uzupełnienia Schura i kiedy jest ona wydajna obliczeniowo, aby z niej korzystać!
Paweł
6

Metoda redukcji modelu

PPR(P)Ax=bkk

P

P=[V][diag(1k)000][WT].

Macierze zasłonięte przez gwiazdy mają znaczenie dla innych rzeczy (takich jak błąd szacowania itp.), Ale na razie unikniemy zajmowania się obcymi szczegółami. Wynika, że

P=VWT

P

Zasadniczo rozwiążesz system

PAx=Pb

VWWTV=IPAx=PbWTy=Vx^x

WTAx^=WTb.

x^Vyx

Dlaczego podejście uzupełniające Schur jest prawdopodobnie lepsze

PAx=bR(P)y=xyyxPxyPAPkAkAA2

VWP. nadal jest „dobrą” matrycą projekcyjną dla wszystkich tych systemów, te dodatkowe koszty prawdopodobnie spowodują, że rozwiązanie zredukowanego systemu będzie droższe niż rozwiązanie oryginalnego systemu.

Wady są bardzo podobne do podejścia JackPoulsona, z tym wyjątkiem, że nie do końca wykorzystujesz wspomnianą strukturę.

Geoff Oxberry
źródło
4

Długa odpowiedź jest ... w pewnym sensie.

Możesz zmienić układ równań tak, aby był jak najdalszy k kolumny to zmienne, dla których chcesz rozwiązać.

Krok 1: Wykonaj eliminację Gaussa, aby matryca była górna trójkątna. Krok 2: rozwiązuj przez podstawienie wsteczne tylko dla pierwszego (ostatniego)k zmienne, którymi jesteś zainteresowany

Pozwoli ci to zaoszczędzić obliczeniową złożoność rozwiązywania ostatnich problemów n-k zmienne poprzez podstawienie wsteczne, które może być tego warte, jeśli njest tak duży, jak mówisz. Należy pamiętać, że w kroku 1 nadal trzeba będzie wykonać sporo pracy.

Pamiętaj również, że ograniczenie kolejności, w której zamierzasz wykonać podstawienie wsteczne, może ograniczyć formę matrycy (zabiera to możliwość wymiany kolumn), co może prowadzić do źle uwarunkowanego systemu, ale ja nie jestem pewny tego - tylko o czym należy pamiętać.

drjrm3
źródło
Wymaga eliminacji gaussowskiej O(n3)) działa, ale podstawienie wstecz wymaga tylko O(n2)). Tak więc, jakn rośnie, procent czasu spędzonego na rozwiązywaniu trójkąta staje się znikomo mały.
Jack Poulson,
dlatego odpowiedź brzmi „w pewnym sensie” zamiast „tak” =)
drjrm3
Ma to sens, że można to zrobić w ten sposób ... Jednak większość obliczeń w eliminacji Gaussa znajduje się w fazie eliminacji do przodu, co powoduje złożoność O (n ^ 3) pomimo skróconej fazy substytucji wstecznej. Miałem nadzieję, że istnieje szybsza metoda ...
Paweł