rozkład LU

W analizie numerycznej i algebrze liniowej , rozkład dolny-górny ( LU ) lub faktoryzacja macierzy jako iloczyn dolnej macierzy trójkątnej i górnej macierzy trójkątnej (patrz rozkład macierzy ) . Produkt czasami zawiera również macierz permutacji . Dekompozycję LU można postrzegać jako postać macierzową eliminacji Gaussa . Komputery zwykle rozwiązują kwadratowe układy równań liniowych za pomocą dekompozycji LU, a także jest kluczowym krokiem podczas odwracania macierzy lub obliczania wyznacznika macierzy . Rozkład LU został wprowadzony przez polskiego matematyka Tadeusza Banachiewicza w 1938 roku. Nazywa się go również rozkładem LR (rozkłada na lewe i prawe trójkątne macierze).

Definicje

Dekompozycja LDU macierzy Walsha

Niech A będzie macierzą kwadratową. Faktoryzacja LU odnosi się do faktoryzacji A , z odpowiednim uporządkowaniem wierszy i / lub kolumn lub permutacji, na dwa czynniki - dolną macierz trójkątną L i górną macierz trójkątną U :

${\ displaystyle A = LU.}$

W dolnej trójkątnej macierzy wszystkie elementy powyżej przekątnej są zerowe, w górnej trójkątnej macierzy wszystkie elementy poniżej przekątnej są zerowe. Na przykład dla macierzy A 3 × 3 jej rozkład LU wygląda następująco:

${\ Displaystyle {\ rozpocząć {bmatrix} a_ {11} i a_ {12} i a_ {13} \\ a_ {21} i a_ {22} i a_ {23} \\ a_ {31} i a_ {32} i a_ {33} \ end{bmatrix}}={\begin{bmatrix}\ell _{11}&0&0\\\ell _{21}&\ell _{22}&0\\\ell _{31}&\ell _{32} &\ell _{33}\end{bmatrix}}{\begin{bmatrix}u_{11}&u_{12}&u_{13}\\0&u_{22}&u_{23}\\0&0&u_{33}\end{ bmacierz}}.}$

Bez odpowiedniego uporządkowania lub permutacji w macierzy faktoryzacja może się nie urzeczywistnić. Na przykład łatwo jest zweryfikować (rozwijając mnożenie macierzy ), że za ${\textstyle a_ {11}=\ell _ {11}u_ {11}$ . Jeśli ${\textstyle a_{11}=0}$ , to co najmniej jeden z ${\textstyle \ ell _ {11}}$ i ${\textstyle u_ {11}}$ musi wynosić zero, co implikuje że albo Ł lub U jest liczbą pojedynczą . Jest to niemożliwe, jeśli A nie jest osobliwe (odwracalne). To jest problem proceduralny. Można go usunąć, po prostu zmieniając kolejność wierszy A , tak aby pierwszy element permutowanej macierzy był różny od zera. Ten sam problem w kolejnych krokach faktoryzacji można usunąć w ten sam sposób; patrz podstawowa procedura poniżej.

Faktoryzacja LU z częściowym obracaniem

Okazuje się, że odpowiednia permutacja w wierszach (lub kolumnach) jest wystarczająca do faktoryzacji LU. Faktoryzacja LU z częściowym obracaniem (LUP) często odnosi się do faktoryzacji LU tylko z permutacjami wierszy:

${\ Displaystyle PA = LU,}$

gdzie L i U są ponownie dolnymi i górnymi macierzami trójkątnymi, a P jest macierzą permutacji , która po pomnożeniu w lewo do A zmienia kolejność wierszy A . Okazuje się, że wszystkie macierze kwadratowe można rozłożyć na czynniki w tej postaci, a faktoryzacja jest w praktyce numerycznie stabilna. To sprawia, że ​​dekompozycja LUP jest użyteczną techniką w praktyce.

Faktoryzacja LU z pełnym obracaniem

Faktoryzacja LU z pełnym obracaniem obejmuje zarówno permutacje wierszy, jak i kolumn:

${\ Displaystyle PAQ = LU,}$

gdzie L , U i P są zdefiniowane jak poprzednio, a Q jest macierzą permutacji, która zmienia kolejność kolumn A .

Dekompozycja dolna-przekątna-górna (LDU).

dolnej przekątnej-górnej (LDU) jest dekompozycją formy

${\ Displaystyle A = LDU,}$

gdzie D jest macierzą diagonalną , a L i U są macierzami jednostkowymi trójkątnymi , co oznacza, że ​​wszystkie wpisy na przekątnych L i U są równe.

Macierze prostokątne

Powyżej wymagaliśmy, aby A była macierzą kwadratową, ale wszystkie te dekompozycje można również uogólnić na macierze prostokątne. W takim przypadku L i D są macierzami kwadratowymi, z których obie mają taką samą liczbę wierszy jak A , a U ma dokładnie takie same wymiary jak A . Górny trójkąt należy interpretować jako mający tylko zero wpisów poniżej głównej przekątnej, która zaczyna się w lewym górnym rogu. Podobnie, bardziej precyzyjnym określeniem dla U jest to, że jest to „forma schodkowa wiersza” macierzy A. _

Przykład

Rozkładamy następującą macierz 2 na 2:

${\ Displaystyle {\ rozpocząć {bmatrix} 4 i 3 \\ 6 i 3 \ koniec {bmatrix}} = {\ rozpocząć {bmatrix} \ ell _ {11} i 0 \\\ ell _ {21} i \ ell _ {22} \ koniec {bmatrix}}{\begin{bmatrix}u_{11}&u_{12}\\0&u_{22}\end{bmatrix}}.}$

Jednym ze sposobów znalezienia rozkładu LU tej prostej macierzy byłoby po prostu rozwiązanie równań liniowych przez sprawdzenie. Rozszerzenie mnożenia macierzy daje

${\ Displaystyle {\ rozpocząć {wyrównane} \ ell _ {11} \ cdot u_ {11} + 0 \ cdot 0 & = 4 \\\ ell _ {11} \ cdot u_ {12} + 0 \ cdot u_ {22} &=3\\\ell _{21}\cdot u_{11}+\ell _{22}\cdot 0&=6\\\ell _{21}\cdot u_{12}+\ell _{22} \cdot u_{22}&=3.\end{wyrównane}}}$

Ten układ równań jest niedookreślony . W tym przypadku dowolne dwa niezerowe elementy L i U są parametrami rozwiązania i można je dowolnie ustawić na dowolną niezerową wartość. Dlatego, aby znaleźć unikalny rozkład LU, konieczne jest nałożenie pewnych ograniczeń na macierze L i U. Na przykład, możemy wygodnie wymagać, aby dolna trójkątna macierz L była jednostkową macierzą trójkątną (tzn. ustawić wszystkie wpisy jej głównej przekątnej na jedynki). Wówczas układ równań ma następujące rozwiązanie:

${\ Displaystyle {\ rozpocząć {wyrównane} \ ell _ {11} = \ ell _ {22} i = 1 \ \\ell _{21}&=1,5\\u_{11}&=4\\u_{12}&=3\\u_{22}&=-1,5\end{wyrównane}}}$

Podstawienie tych wartości do powyższego rozkładu LU daje wyniki

${\ Displaystyle {\ rozpocząć {bmatrix} 4 i 3 \\ 6 i 3 \ koniec {bmatrix}} = {\ rozpocząć {bmatrix} 1 i 0 \\ 1,5 i 1 \ koniec {bmatrix}} {\ rozpocząć {bmatrix} 4 i 3 \\ 0 i -1,5 \ koniec {bmacierz}}.}$

Istnienie i wyjątkowość

Macierze kwadratowe

$Dowolna$ macierz kwadratowa rozkłady na czynniki LUP i PLU . Jeśli $odwracalny$ , to dopuszcza rozkład na czynniki LU (lub LDU ) wtedy i tylko wtedy, gdy wszystkie jego wiodące główne drugorzędne są niezerowe (na przykład [ ${\ Displaystyle {\ rozpocząć {bmatrix} 0 i 1$ nie dopuszcza LU ani LDU faktoryzacja). Jeśli ${\ textstyle k}$ jest macierzą liczby pojedynczej o randze $k}$$\ textstyle$ , to dopuszcza faktoryzację LU , jeśli pierwsze główne drugorzędne są niezerowe, chociaż odwrotność nie jest prawdziwa.

Jeśli kwadratowa, odwracalna macierz ma LDU (rozkład na czynniki ze wszystkimi przekątnymi L i U równymi 1), to faktoryzacja jest niepowtarzalna. W takim przypadku LU jest również wyjątkowa $)$ jeśli wymagamy, aby przekątna ${\textstyle L}$ (lub składała się z jedynek.

Ogólnie rzecz biorąc, każda macierz kwadratowa może mieć jedną z następujących cech: ZA $A_ {n \ razy n}}$

1. unikalna faktoryzacja LU (jak wspomniano powyżej);
2. nieskończenie wiele faktoryzacji LU, jeśli dwie lub więcej dowolnych pierwszych ( n -1) kolumn jest liniowo zależnych lub dowolna z pierwszych ( n -1) kolumn ma wartość 0;
3. brak faktoryzacji LU, jeśli pierwsze ( n −1) kolumny są niezerowe i liniowo niezależne, a co najmniej jedna wiodąca główna podrzędna ma wartość zero.

W przypadku 3 można przybliżyć faktoryzację LU, zmieniając wpis ukośny na za $Displaystyle a_ {jj}}$ na ${\ Displaystyle a_ {jj} \ pm \ varepsilon},$ aby uniknąć głównego zera na początku drobny.

Symetryczne macierze dodatnio określone

Jeśli A jest symetryczną (lub hermitowską , jeśli A jest zespolona) dodatnio określoną macierzą , możemy tak ułożyć sprawy, że U jest sprzężoną transpozycją L. Oznacza to, że możemy zapisać A jako

${\ Displaystyle A = LL ^ {*}. \,}$

Ten rozkład nazywa się rozkładem Choleskiego . Rozkład Cholesky'ego zawsze istnieje i jest unikalny — pod warunkiem, że macierz jest dodatnio określona. Ponadto obliczenie dekompozycji Cholesky'ego jest bardziej wydajne i liczbowo bardziej stabilne niż obliczenie niektórych innych dekompozycji LU.

Macierze ogólne

W przypadku macierzy (niekoniecznie odwracalnej) w dowolnym polu znane są dokładne warunki konieczne i wystarczające, w których ma ona rozkład na czynniki LU. Warunki są wyrażone jako rangi pewnych podmacierzy. Algorytm eliminacji Gaussa służący do uzyskiwania rozkładu LU został również rozszerzony na ten najbardziej ogólny przypadek.

Algorytmy

Zamknięta formuła

Gdy faktoryzacja LDU istnieje i jest niepowtarzalna, istnieje zamknięty (jawny) wzór na elementy L , D i U pod względem stosunków wyznaczników pewnych podmacierzy oryginalnej macierzy A . W szczególności ${\textstyle D_ {1}=A_ {1,1}}$ i dla ${\textstyle i=2,\ldots, n}$ , ${\ textstyle D_ {i}}$$_$ stosunkiem $podmacierzy$ do głównej Obliczanie wyznaczników jest kosztowne obliczeniowo , więc ten jawny wzór nie jest używany w praktyce.

Korzystanie z eliminacji Gaussa

Poniższy algorytm jest zasadniczo zmodyfikowaną formą eliminacji Gaussa . $Obliczenie$ rozkładu LU przy użyciu tego algorytmu wymaga niższego rzędu. Częściowe obracanie dodaje tylko składnik kwadratowy; nie dotyczy to pełnego obrotu.

Procedura

Biorąc pod uwagę macierz N × N $\ Displaystyle A = (a_ {i, j}) _ {1 \ równoważnik ja, j \ równoważnik N}}$ , zdefiniuj $}$$}$ jako macierz , w której niezbędne wiersze zostały zamienione, aby spełnić pożądane warunki (takie jak częściowe obrócenie) dla pierwszej kolumny. Indeks górny w nawiasie (np. ${\ Displaystyle (0)}$ ) macierzy jest wersją macierzy ${\ Displaystyle A} .$ Macierz jest $\$ , w której elementy poniżej głównej przekątnej zostały $) {$ wyeliminowane do 0 poprzez eliminację Gaussa dla pierwszego $n}$$Displaystyle A ^ {($ kolumn, a niezbędne wiersze zostały zamienione, aby spełnić żądane warunki dla ${\ Displaystyle (n + 1) ^ {th}}$ kolumna.

Wykonujemy operację $}}$${\ Displaystyle row_ {i} = row_ {i} - (\ ell _ {i, n}) \ cdot row_ {$$z$ dla każdego wiersza elementami (oznaczonymi jako ${\ Displaystyle a_ {i, n} ^ {(n-1)}}$ gdzie ${\ Displaystyle i = n + 1, \ kropka, N}$ ) poniżej głównej przekątnej w n -tej kolumnie ${\ Displaystyle A ^ {(n-1) }}$ . Do tej operacji

${\ Displaystyle \ ell _ {i, n}: = {\ Frac {a_ {i, n} ^ {(n-1)}} {a_ {n, n} ^ {(n-1)}}}. }$

Wykonujemy te operacje na wierszach, aby wyeliminować elementy do zera. ${\ Displaystyle a_ {i, n} ^ {(n-1)}}$$Po odjęciu$ możemy zamienić wiersze, aby dla Możemy tutaj zamienić wiersze, aby wykonać częściowe przestawienie lub ponieważ element ${\ Displaystyle a_ {n + 1, n + 1}}$ na głównej przekątnej wynosi zero (a zatem nie można go użyć do implementacji eliminacji Gaussa).

Definiujemy ostateczną macierz permutacji $co$ macierz tożsamości, która ma wszystkie te same wiersze zamienione w tej samej kolejności $,$ macierz

Po przeprowadzeniu operacji na wierszach dla pierwszych kolumn otrzymaliśmy górną macierz trójkątną $N-1)}}$ która wynosi $A$ oznaczony przez ${\ Displaystyle U}$ . Uwaga, możemy oznaczyć jako jako $\ Displaystyle$ tą kolumnę $A ^ {(N-1)}}$${\ Displaystyle A ^ {(N-1)}}$ nie ma warunków, dla których należy zamienić wiersze. Możemy również obliczyć $ell$$oznaczoną$ trójkątną $jako$ , tak że , bezpośrednio wprowadzając wartości wartości za pomocą poniższego wzoru.

${\ Displaystyle L = {\ rozpocząć {pmatrix} 1&&&&&\\\ ell _ {2,1}&\ ddots &&&&\\&\ ddots & 1&&&\\\ vdots &\ cdots &\ ell _ {n + 1, n} &1&&\\&&\vdots &\ddots &\ddots &\\\ell _{N,1}&\cdots &\ell _{N,n}&\cdots &\ell _{N,N-1}&1 \end{pmacierz}}}$

Przykład

Jeśli otrzymamy macierz

$${\ Displaystyle A = {\ rozpocząć {pmatrix} 0 i 5 i {\ Frac {22} {3}} \\ 4 i 2 i 1 \\ 2 i 7 i 9 \\\ koniec {pmatrix}} ,}$$

zdecydujemy się zaimplementować częściowe obracanie, a tym samym zamienić pierwszy i drugi wiersz, tak aby nasza macierz $i$ iteracja naszej macierzy stały się odpowiednio ${\ displaystyle P}$

$${\ Displaystyle A ^ {(0)} = {\ początek {pmatrix} 4 i 2 i 1 \\ 0 i 5 & {\ Frac {22} {3}} \\ 2 i 7 i 9 \\\ koniec {pmatrix}}, \ quad P ^ {(0 )}={\początek{pmacierzy}0&1&0\\1&0&0\\0&0&1\\\koniec{pmacierzy}}.}$$

Po zamianie rzędów możemy wyeliminować elementy znajdujące się poniżej głównej przekątnej w pierwszej kolumnie, wykonując czynności

$${\ Displaystyle {\ rozpocząć {wyrównane }{0}wiersz_{2}=wiersz_{2}-(\ell _{2,1})\cdot row_{1}\\row_{3}=wiersz_{3}-(\ell _{3,1 })\cdot wiersz_{1}\end{wyrównanydat}}}$$

takie, że

$${\ Displaystyle {\ rozpocząć {wyrównane} {0} \ ell _ {2,1} = {\ Frac {0} {4}} = 0 \ \\ell _{3,1}={\frac {2}{4}}=0.5\end{wyrównanydat}}}$$

Po odjęciu tych wierszy wyprowadziliśmy z macierzy ${\ Displaystyle A ^ {(0)}}$

$${\ Displaystyle {\ rozpocząć {pmatrix} 4 i 2 i 1 \\ 0 i 5 i {\ Frac {22} {3}} \\ 0 i 6 i 8,5 \\\ koniec {pmatrix}}.}$$

Ponieważ wdrażamy częściowe obracanie, zamieniamy odpowiednio drugi i trzeci wiersz naszej macierzy pochodnej oraz bieżącą wersję naszej macierzy, aby uzyskać ${\ displaystyle P}$

$${\ Displaystyle A ^ {(1)} = {\ początek {pmatrix} 4 i 2 i 1 \\ 0 i 6 i 8,5 \\ 0 i 5 i {\ Frac {22} {3}} \\\ koniec {pmatrix}}, \ quad P ^ { (1)}={\początek{pmacierzy}0&1&0\\0&0&1\\1&0&0\\\koniec{pmacierzy}}.}$$

Teraz eliminujemy elementy poniżej głównej przekątnej w drugiej kolumnie, wykonując ${\ displaystyle row_ {3} = row_ {3} - (\ ell _ {3,2}) \ cdot row_ {2}}$ takie, że ${\ textstyle \ ell _ {3,2} = {\ Frac {5} {6}}}$ . Ponieważ żadne elementy niezerowe nie istnieją poniżej głównej przekątnej w naszej bieżącej iteracji ${\ displaystyle$$wierszy$ , to odejmowanie wierszy wyprowadza naszą ostateczną macierz (oznaczoną jako ) $A$ } i końcowa macierz: ${\ displaystyle P}$

$${\ Displaystyle A ^ {(2)} = A ^ {(N-1)} = U = {\ początek {pmatrix} 4 i 2 i 1 \\ 0 i 6 i 8,5 \\ 0 i 0 i 0,25 \\\ koniec {pmatrix}}, \ quad P={\początek{pmacierz}0&1&0\\0&0&1\\1&0&0\\\koniec{pmacierz}}.}$$

Teraz możemy uzyskać naszą ostateczną macierz: ${\ displaystyle L}$

$${\ Displaystyle L = {\ rozpocząć {pmatrix} 1 i 0 i 0 \\\ ell _ {3,1} &1&0\\\ell _{2,1}&\ell _{3,2}&1\\\end{pmatrix}}={\begin{pmatrix}1&0&0\\0.5&1&0\\0&{\frac {5} {6}}&1\\\koniec{pmacierzy}}}$$

Teraz te macierze mają taką relację, że ${\ displaystyle PA = LU}$ .

Relacje, gdy żadne wiersze nie są zamienione

Jeśli w ogóle nie zamieniliśmy wierszy podczas tego procesu, możemy wykonać operacje na wierszach jednocześnie dla każdej kolumny, $L$$^ {(n-1)},}$ ) $A ^$ { (n)}:= L_ { n} $N z$ gdzie jest macierzą identyczności N × n -tą kolumną zastąpioną wektor ${\ Displaystyle {\ rozpocząć {pmatrix} 0 & \ dotsm & 0 i 1 & - \ ell _ {n + 1, n} & \ kropka & - \ ell _ {N, n} \ koniec {pmatrix}} ^ {\ textsf {T} }.}$ Innymi słowy, dolna macierz trójkątna

${\ Displaystyle L_ {n} = {\ rozpocząć {pmatrix} 1&&&&&\\&\ ddots &&&&\\&&1&&&\\&&-\ ell _ {n+1,n}&&&\\&&\vdots &&\ddots &\\ &&-\ell _{N,n}&&&1\end{pmacierz}}.}$

$\ Displaystyle A ^ {(n)}: = L_$ na wierszach dla pierwszych kolumn za pomocą $1$ jest równoważna znalezieniu rozkładu

$${\ Displaystyle A = L_ {1} ^ {-1} L_ {1} A ^ {(0)} = L_ {1} ^ {- 1} A ^ {(1)} = L_ {1} ^ {- 1}L_{2}^{-1}L_{2}A^{(1)}=L_{1}^{-1}L_{2}^{-1}A^{(2)}=\ dotsm =L_{1}^{-1}\dotsm L_{N-1}^{-1}A^{(N-1)}.}$$

Oznacz ${\ textstyle L = L_ {1} ^ {- 1} \ kropka L_ {N-1} ^ {-1}}$ tak, że ${\ Displaystyle A = LA ^ {(N-1)} = LU}$ .

Teraz obliczmy sekwencję ${\ Displaystyle L_ {1} ^ {-1} \ kropka L_ {N-1} ^ {-1}}$ . Wiemy $ma$ następującą

${\ Displaystyle L_ {n} ^ {- 1} = {\ rozpocząć {pmatrix} 1 &&&&&\\&\ ddots &&&&\ \&&1&&&\\&&\ell _{n+1,n}&&&\\&&\vdots &&\ddots &\\&&\ell _{N,n}&&&1\end{pmatrix}}}$

Jeśli istnieją dwie dolne trójkątne macierze z jedynkami na głównej przekątnej i żadna z nich nie ma elementu niezerowego poniżej głównej przekątnej w tej samej kolumnie co druga, możemy uwzględnić wszystkie elementy niezerowe w tym samym miejscu w produkcie z dwóch macierzy. Na przykład:

${\ Displaystyle \ lewo ({\ rozpocząć {tablica} {ccccc} 1&0&0&0&0 \\77&1&0&0&0 \\12&0&1&0&0\\63&0&0&1&0\\7&0&0&0&1\koniec {tablica}}\prawo)\lewo ({\rozpocząć {tablica}}1&0&0&0&0 \\0&1&0&0&0\\0&22&1&0&0\\0&33&0&1&0\\0&44&0&0&1\end{array}}\right)=\left({\begin{array}{ccccc}1&0&0&0&0\\77&1&0&0&0\\12&22&1&0&0\\63&33&0&1&0\ \7&44&0&0&1\koniec{ tablica}}\prawo)}$

$wcześniej$ pomnóż razem i skondensowaną macierz oznaczoną jako jak wspomniano $.$ Używając macierzy ${\textstyle L}$ , otrzymujemy ${\ Displaystyle A = LU.}$

Oczywiste jest, że aby ten algorytm działał, na każdym kroku trzeba mieć ${\ Displaystyle a_ {n, n} ^ {(n-1)} \ neq 0}$ (zobacz definicję ${\ Displaystyle \ ell _ {i, n}}$ ). Jeśli to założenie zawiedzie w pewnym momencie, przed kontynuowaniem należy zamienić n -ty wiersz z innym wierszem poniżej. To dlatego rozkład LU ogólnie wygląda tak: ${\ Displaystyle P ^ {- 1} A = LU}$ .

LU Rozkład krouta

Zauważ, że rozkład uzyskany dzięki tej procedurze jest rozkładem Doolittle'a : główna przekątna L składa się wyłącznie z 1 s. Gdybyśmy postępowali usuwając elementy powyżej głównej przekątnej, dodając wielokrotności kolumn ( zamiast usuwać elementy poniżej przekątnej, dodając wielokrotności wierszy ) , otrzymalibyśmy rozkład Crouta , w którym główna przekątna U wynosi 1 s .

Innym (równoważnym) sposobem uzyskania rozkładu Crouta danej macierzy A jest uzyskanie rozkładu Doolittle'a transpozycji A . $sekcji$ , jeśli rozkładem LU uzyskanym za pomocą algorytmu przedstawionego w tej ${\textstyle L=U_ {0}^{\textsf {T}}}$ i ${\textstyle U=L_{0}^{\textsf {T}}}$ mamy to ${\ displaystyle A = LU}$ jest rozkładem Crouta.

Poprzez rekurencję

Cormen i in. opisz rekurencyjny algorytm dekompozycji LUP.

Mając macierz A , niech P ₁ będzie macierzą permutacji taką, że

${\ Displaystyle P_ {1} A = \ lewo ({\ rozpocząć {tablica} {c | ccc} a & & w ^ {\ textsf {T}} & \\\ hline &&& \\v&&A'&\\&&&\end{tablica}}\right)}$ ,

gdzie ${\ textstyle a\ neq 0}$ , jeśli w pierwszej kolumnie A znajduje się wpis niezerowy ; lub w przeciwnym razie weź P ₁ jako macierz tożsamości. Niech teraz ${\textstyle c=1/a}$ , jeśli ${\textstyle a\neq 0}$ ; lub ${\ textstyle c=0}$ w przeciwnym razie. Mamy

${\ Displaystyle P_ {1} A = \ lewo ({\ początek {tablica} {c | ccc} 1&& 0& \\\ hline &&&\\cv&&I_ {n-1}&\\&&&\ koniec {tablica}} \ prawej) \left({\begin{array}{c|c}a&w^{\textsf {T}}\\\hline &\\0&A'-cvw^{\textsf {T}}\\&\end{array} }\Prawidłowy).}$

$_ U'}$ możemy . Niech ${\ textstyle v'= P'v}$ . Dlatego

${\ Displaystyle \ lewo ({\ początek {tablica} {c|ccc} 1&&0&\\\hline &&&\\0&&P'&\\&&&\koniec {tablica}}\right)P_{1}A=\left({{ \begin{array}{c|ccc}1&&0&\\\hline &&&\\cv'&&L'&\\&&&\end{array}}\right)\left({\begin{array}{c|ccc}a&&w ^{\textsf {T}}&\\\hline &&&\\0&&U'&\\&&&\end{tablica}}\right),}$

który jest rozkładem LUP A .

Algorytm losowy

Możliwe jest znalezienie przybliżenia niskiego rzędu do rozkładu LU przy użyciu losowego algorytmu. Biorąc pod uwagę macierz wejściową $P$ pożądaną niską rangę $\$ losowa jednostka logiczna zwraca macierze permutacji $textstyle P, Q}$ i dolne/górne macierze trapezowe ${\ textstyle L, U}$ o rozmiarze ${\ textstyle m \ razy k}$ i ${\ textstyle k \ razy n}$ odpowiednio, takie, że z dużym prawdopodobieństwem ${\textstyle \left\|PAQ-LU\right\|_{2}\leq C\sigma _{k+1 }}$ , gdzie ${\ textstyle (k +1)}$ stałą zależną od parametrów algorytmu i $\ textstyle \ sigma _ {k + 1}}$ to $+$ -ta pojedyncza wartość macierzy wejściowej ${\textstyle A}$ .

Złożoność teoretyczna

Jeśli dwie macierze rzędu n można pomnożyć w czasie M ( n ), gdzie M ( n ) ≥ n ^a dla pewnego a > 2, to rozkład LU można obliczyć w czasie O ( M ( n )). Oznacza to na przykład, że istnieje algorytm O( n ^2,376 ) oparty na algorytmie Coppersmitha-Winograda .

Dekompozycja macierzy rzadkiej

Opracowano specjalne algorytmy do rozkładu na czynniki dużych macierzy rzadkich . Algorytmy te próbują znaleźć rzadkie czynniki L i U . W idealnym przypadku koszt obliczeń jest określany przez liczbę niezerowych wpisów, a nie przez rozmiar macierzy.

Algorytmy te wykorzystują swobodę wymiany wierszy i kolumn, aby zminimalizować wypełnianie (wpisy, które zmieniają się od początkowego zera do wartości niezerowej podczas wykonywania algorytmu).

Ogólne traktowanie porządków, które minimalizują wypełnienie, można rozwiązać za pomocą teorii grafów .

Aplikacje

Rozwiązywanie równań liniowych

Dany układ równań liniowych w postaci macierzowej

${\ Displaystyle A \ mathbf {x} = \ mathbf {b},}$

chcemy rozwiązać równanie dla x , biorąc pod uwagę A i b . Załóżmy ${\$$textstyle LU\mathbf {x} = P\mathbf {$ uzyskaliśmy już rozkład LUP A taki, że , więc .

W tym przypadku rozwiązanie odbywa się w dwóch logicznych krokach:

1. Najpierw rozwiązujemy równanie ${\textstyle L\mathbf {y} = P\mathbf {b}}$ dla y .
2. Po drugie, rozwiązujemy równanie dla x ${\textstyle U\mathbf {x} =\mathbf {y}$ }

W obu przypadkach mamy do czynienia z macierzami trójkątnymi ( L i U ), które można rozwiązać bezpośrednio przez podstawienie w przód iw tył bez stosowania procesu eliminacji Gaussa (jednak potrzebujemy tego procesu lub jego odpowiednika do obliczenia samego rozkładu LU ).

Powyższą procedurę można wielokrotnie zastosować do rozwiązania równania wiele razy dla różnych b . W takim przypadku szybsze (i wygodniejsze) jest wykonanie dekompozycji LU macierzy A raz, a następnie rozwiązanie trójkątnych macierzy dla różnych b , zamiast stosowania eliminacji Gaussa za każdym razem. Można by pomyśleć, że macierze L i U „zakodowały” proces eliminacji Gaussa.

Koszt rozwiązania układu równań liniowych $wynosi$$jeśli$ macierz ma ${\ textstyle n}$ . To sprawia $,$ dwa razy szybszy niż algorytmy oparte na dekompozycji QR co kosztuje około operacji gdy Householder odzwierciedla są używane. Z tego powodu zwykle preferowana jest dekompozycja LU.

Odwracanie macierzy

Przy rozwiązywaniu układów równań b jest zwykle traktowane jako wektor o długości równej wysokości macierzy A . Jednak w odwróceniu macierzy zamiast wektora b mamy macierz B , gdzie B jest macierzą n -na- p , więc próbujemy znaleźć macierz X (również macierz n -na- p ):

${\ Displaystyle AX = LUX = B.}$

Możemy użyć tego samego algorytmu przedstawionego wcześniej do rozwiązania dla każdej kolumny macierzy X . Załóżmy teraz, że B jest macierzą tożsamości o rozmiarze n . Wynikałoby z tego, że wynik X musi być odwrotnością A .

Obliczanie wyznacznika

Biorąc $A$ rozkład LUP wyznacznik A można bezpośrednio

${\ Displaystyle \ det (A) = \ det \ lewo (P ^ {-1} \ prawo) \ det (L) \ det (U) = (-1) ^ {S} \ lewo (\ prod _ {i =1}^{n}l_{ii}\right)\left(\prod _{i=1}^{n}u_{ii}\right).}$

Drugie równanie wynika z faktu, że wyznacznik macierzy trójkątnej jest po prostu iloczynem jej przekątnych wpisów, a wyznacznik macierzy permutacji jest równy (−1) ^S , gdzie S jest liczbą wymian wierszy w rozkładzie .

W przypadku rozkładu LU z pełnym obrotem, równa się również prawej stronie powyższego równania, jeśli pozwolimy $całkowitą$ być wymian wierszy i

Ta sama metoda może być łatwo zastosowana do rozkładu jednostek logicznych przez ustawienie P równego macierzy tożsamości.

Przykłady kodu

Przykład kodu C

 /* WEJŚCIE: A - tablica wskaźników do wierszy macierzy kwadratowej o wymiarze N  * Tol - mała tolerancja do wykrycia awarii, gdy macierz jest bliska degeneracji  * WYJŚCIE: Macierz A jest zmieniona, zawiera kopię obu macierzy LE i U jako A=(LE)+U takie, że P*A=L*U.  * Macierz permutacji nie jest przechowywana jako macierz, ale w wektorze całkowitym P o rozmiarze N+1  * zawierającym indeksy kolumn, gdzie macierz permutacji ma „1”. Ostatni element P[N]=S+N,   * gdzie S to liczba wymian wierszy potrzebnych do obliczenia wyznacznika, det(P)=(-1)^S  */  int          

         
       

       0    
         LUPDecompose  (  podwójne  **  A  ,  int  N  ,  podwójne  Tol  ,  int  *  P  )  {  int  i  ,  j  ,  k  ,  imax  ;  podwójne  maxA  ,  *  ptr  ,  absA  ;  dla  (  ja  =  ;  ja  <=  N  ;  ja  ++  )  P  [  ja    

       0     
          
          

               
                ]  =  ja  ;  // Macierz permutacji jednostek, P[N] zainicjowana N  for  (  i  =  ;  i  <  N  ;  i  ++  )  {  maxA  =  0.0  ;  imax  =  ja  ;  for  (  k  =  i  ;  k  <  N  ;  k  ++  )  if  ((  absA  =  fabs  (  A  [     
                  
                  
            

             0 

            
            
              
               k  ][  ja  ]))  >  maxA  )  {  maxA  =  absA  ;  imaks  =  k  ;  }  jeśli  (  maxA  <  Tol  )  powrót  ;  //niepowodzenie, macierz jest zdegenerowana  if  (  imax  !=  i  )  {  //obrót P  j  =  P  [  i  ];  P.  [  ja  ]  =  P.  [ 
              

            
              
              
              

            
            
        

             imax  ];  P.  [  imax  ]  =  j  ;  //przestawne wiersze A  ptr  =  A  [  i  ];  ZA  [  ja  ]  =  ZA  [  imaks  ];  ZA  [  imax  ]  =  ptr  ;  //liczenie osi zaczynając od N (dla wyznacznika)  P  [  N  ]  ++  ;  }  dla  (  j  =  i  +       
              

                     
                    1  ;  j  <  N  ;  j  ++  )  {  ZA  [  j  ][  ja  ]  /=  ZA  [  ja  ][  ja  ];  dla  (  k  =  ja  +  1  ;  k  <  N  ;  k  ++  )  ZA  [  jot  ][  k  ]  -=  ZA  [  jot  ][  ja  ]  *  
        
    

       





           

      ZA  [  ja  ][  k  ];  }  }  powrót  1  ;  //dekompozycja zakończona  }  /* INPUT: A,P wypełnione w LUPDecompose; b - wektor prawy;  N - wymiar   * WYJŚCIE: x - wektor rozwiązania A*x=b  */  void  LUPSolve  (  double  **  A  ,  int  *  P  ,  double  *  b  ,  int  N  ,  double  *  x  )  {  for  (    0     
          

            0    
                 int  ja  =  ;  ja  <  N  ;  ja  ++  )  {  x  [  ja  ]  =  b  [  P  [  ja  ]];  dla  (  int  k  =  ;  k  <  ja  ;  k  ++  )  x  [  ja  ]  -=  ZA  [  ja  ][  k  ]  *  x  [  k  ]; 
    

             0  
                  
                

         }  for  (  int  ja  =  N  -  1  ;  ja  >=  ;  ja  --  )  {  for  (  int  k  =  ja  +  1  ;  k  <  N  ;  k  ++  )  x  [  ja  ]  -=  ZA  [  ja  ][  k  ]  *  x  [  k  ];  x  [  ja   
    





         
  
        0     ]  /=  A  [  ja  ][  ja  ];  }  }  /* INPUT: A,P wypełnione w LUPDecompose; N - wymiar   * WYJŚCIE: IA jest odwrotnością początkowej macierzy  */  void  LUPInvert  (  double  **  A  ,  int  *  P  ,  int  N  ,  double  **  IA  )  {  for  (  int  j  =  ;  j  <  N  ;  j  
            0     
                    

                0    
                 ++  )  {  dla  (  int  ja  =  ;  ja  <  N  ;  ja  ++  )  {  JA  [  ja  ][  jot  ]  =  P.  [  ja  ]  ==  jot  ?  1,0  :  0,0  ;  dla  (  int  k  =  ;  k  <  ja  ;  k  ++  )  IA  [  ja     
        

                 0  
                      
                 ][  jot  ]  -=  ZA  [  ja  ][  k  ]  *  IA  [  k  ][  jot  ];  }  for  (  int  ja  =  N  -  1  ;  ja  >=  ;  ja  --  )  {  for  (  int  k  =  ja  +  1  ;  k  <  N  ;  k  ++  )  IA  [     

              
        
    





   ja  ][  j  ]  -=  ZA  [  ja  ][  k  ]  *  IA  [  k  ][  j  ];  IA  [  ja  ][  j  ]  /=  A  [  ja  ][  ja  ];  }  }  }  /* INPUT: A,P wypełnione w LUPDecompose; N - wymiar.   * WYJŚCIE: Funkcja zwraca wyznacznik macierzy początkowej  */  double  LUPDeterminant  (  double  **      

       00

            
          

           ZA  ,  int  *  P  ,  int  N  )  {  podwójne  det  =  ZA  [  ][  ];  for  (  int  i  =  1  ;  ja  <  N  ;  ja  ++  )  det  *=  ZA  [  ja  ][  ja  ];  powrót  (  P  [  N  ]  -  N  )  %  2  ==  0    
 ?  det  :  -  det  ;  } 

Przykład kodu C#

  

           
    
        
             
           0
            0 public  class  SystemOfLinearEquations  {  public  double  []  SolveUsingLU  (  double  [,]  matrix  ,  double  []  rightPart  ,  int  n  )  {  // dekompozycja macierzy  double  [,]  lu  =  new  double  [  n  ,  n  ];  podwójna  suma  =  ;  dla  (  int  i  =  ;     
        
                    
            
                  0
                    0    
                           ja  <  n  ;  ja  ++)  {  dla  (  int  j  =  ja  ;  j  <  n  ;  j  ++)  {  suma  =  ;  for  (  int  k  =  ;  k  <  ja  ;  k  ++)  suma  +=  lu  [  ja  ,  k  ]  *  lu  [  k  ,  j  ]; 
                      
            
                      
            
                  0
                    0    
                       lu  [  ja  ,  j  ]  =  macierz  [  ja  ,  j  ]  -  suma  ;  }  for  (  int  j  =  i  +  1  ;  j  <  n  ;  j  ++)  {  suma  =  ;  for  (  int  k  =  ;  k  <  i  ;  k  ++)  suma  +=  lu     
                           
            
        

        
        
             [  j  ,  k  ]  *  lu  [  k  ,  ja  ];  lu  [  j  ,  ja  ]  =  (  1  /  lu  [  ja  ,  ja  ])  *  (  macierz  [  j  ,  ja  ]  -  suma  );  }  }  // lu = L+UI  // znajdź rozwiązanie Ly = b  double  []  y  =  new  double 
            0    
        
              0
                0    
                     
               [  n  ];  for  (  int  ja  =  ;  ja  <  n  ;  ja  ++)  {  suma  =  ;  for  (  int  k  =  ;  k  <  ja  ;  k  ++)  suma  +=  lu  [  ja  ,  k  ]  *  y  [  k  ];  y  [  i  ]  =  prawa część  [   
        
        
            
                 0 
        
              0
                       ja  ]  -  suma  ;  }  // znajdź rozwiązanie Ux = y  double  []  x  =  new  double  [  n  ];  for  (  int  ja  =  n  -  1  ;  ja  >=  ;  ja  --)  {  suma  =  ;  dla  (  int  k  =  i  +  1  ;  k  <  n  ;  k  ++) 
                     
                     
        
         
    
 suma  +=  lu  [  ja  ,  k  ]  *  x  [  k  ];  x  [  ja  ]  =  (  1  /  lu  [  ja  ,  ja  ])  *  (  y  [  ja  ]  -  suma  );  }  zwróćx  ;  _  }  } 

Przykład kodu MATLAB

   
      
      
    
       
             
               funkcja  LU  =  LUDecompDoolittle  (  A  )  n  =  długość  (  A  );  LU  =  A  ;  % rozkładu macierzy, Metoda Doolittle'a  dla  i  =  1  :  1  :  n  dla  j  =  1  :(  i  -  1  )  LU  (  i  ,  j  )  =  (  LU  (  i  ,         
        
          
            j  )  -  LU  (  ja  ,  1  :(  j  -  1  ))  *  LU  (  1  :(  j  -  1  ),  j  ))  /  LU  (  j  ,  j  );  koniec  j  =  ja  :  n  ;  LU  (  ja  ,  j  )  =  LU  (  ja  ,  j  )  -      
    
    


   
      
      
     LU  (  ja  ,  1  :(  ja  -  1  ))  *  LU  (  1  :(  ja  -  1  ),  j  );  koniec %LU =   funkcja  końcowa  L+UI  x  =  SolveLinearSystem  (  LU, B  )  n  =  długość  (  LU  );  y  =  zera  (  rozmiar  (  B  ));  % znajdź rozwiązanie Ly = B 
       
            
    
    
      
        dla  i  =  1  :  n  y  (  ja  ,:)  =  B  (  ja  ,:)  -  LU  (  ja  ,  1  :  ja  )  *  y  (  1  :  ja  ,:);  end  % znajdź rozwiązanie Ux = y  x  =  zera  (  size  (  B  ));  dla  i  =  n  :(  -  1  ):  1 
                 
        


           x  (  ja  ,:)  =  (  y  (  ja  ,:)  -  LU  (  ja  , (  ja  +  1  ):  n  )  *  x  ((  ja  +  1  ):  n  ,:))  /  LU  (  ja  ,  ja  );  koniec  koniec  A  =  [  4  3  3  ;  6  3  3  ;  3  4   
  
               
   
   3  ]  LU  =  LUDecomp Doolittle  (  A  )  B  =  [  1  2  3  ;  4  5  6  ;  7  8  9  ;  10  11  12  ]  '  x  =  Rozwiąż układ liniowy  (  LU  ,  B  )  A  *  x 

Zobacz też

• Zablokuj dekompozycję LU
• rozkład Bruhata
• Rozkład Choleskiego
• Rozkład matrycy krouta
• Niepełna faktoryzacja jednostek logicznych
• Redukcja LU
• Dekompozycja macierzy
• rozkład QR

Notatki

•    Bunch, James R.; Hopcroft, John (1974), „Trójkątna faktoryzacja i inwersja przez szybkie mnożenie macierzy”, Mathematics of Computation , 28 (125): 231–236, doi : 10,2307/2005828 , ISSN 0025-5718 , JSTOR 2005828 .
•   Cormen, Thomas H .; Leiserson, Charles E .; Rivest, Ronald L .; Stein, Clifford (2001), Wprowadzenie do algorytmów , MIT Press i McGraw-Hill, ISBN 978-0-262-03293-3 .
•   Golub, Gene H .; Van Loan, Charles F. (1996), Matrix Obliczenia (wyd. 3), Baltimore: Johns Hopkins, ISBN 978-0-8018-5414-9 .
•   Róg, Roger A.; Johnson, Charles R. (1985), Analiza macierzy , Cambridge University Press, ISBN 978-0-521-38632-6 . Patrz sekcja 3.5. N -1
•   Właściciel domu, Alston S. (1975), Teoria macierzy w analizie numerycznej , New York: Dover Publications , MR 0378371 .
• Okuniew, Paweł; Johnson, Charles R. (1997), Niezbędne i wystarczające warunki istnienia LU Faktoryzacja dowolnej macierzy , arXiv : math.NA/0506382 .
•   Poole, David (2006), Linear Algebra: A Modern Introduction (wyd. 2), Kanada: Thomson Brooks / Cole, ISBN 978-0-534-99845-5 .
•   Prasa, WH; Teukolsky, SA; Vetterling, WT; Flannery, BP (2007), „Sekcja 2.3” , Przepisy numeryczne: sztuka obliczeń naukowych (wyd. 3), Nowy Jork: Cambridge University Press, ISBN 978-0-521-88068-8 .
•   Trefethen, Lloyd N .; Bau, David (1997), Numeryczna algebra liniowa , Filadelfia: Towarzystwo Matematyki Przemysłowej i Stosowanej, ISBN 978-0-89871-361-9 .
• Rigotti, Luca (2001), ECON 2001 - Wprowadzenie do metod matematycznych, Wykład 8

Linki zewnętrzne

Bibliografia

• Dekompozycja LU na MathWorld .
• Dekompozycja LU na Math-Linux .
• Dekompozycja LU w Holistycznym Instytucie Metod Numerycznych
• Faktoryzacja macierzy LU . Odniesienie do MATLAB-a.

Kod komputerowy

• LAPACK to zbiór podprogramów FORTRAN do rozwiązywania problemów z gęstą algebrą liniową
• ALGLIB zawiera częściowy port LAPACK do C++, C#, Delphi itp.
• Kod C++ , prof. J. Loomis, University of Dayton
• Kod C , biblioteka źródeł matematycznych
• Kod rdzy
• LU w X10

Zasoby online

• WebApp opisowo rozwiązujący układy równań liniowych za pomocą rozkładu LU
• Kalkulator macierzowy z krokami, w tym dekompozycją LU ,
• Narzędzie dekompozycji LU , uni-bonn.de
• Dekompozycja LU autorstwa Eda Pegga, Jr. , The Wolfram Demonstrations Project , 2007.