Rozkład LU - LU decomposition
W analizie numerycznej i liniowego Algebra , niższych górnej ( LU ) rozkładu lub faktoryzacji współczynniki A macierz jako produkt o niższej macierzy trójkątnej i górnej trójkątnej matrycy. Produkt czasami zawiera również macierz permutacji . Dekompozycja LU może być postrzegana jako macierzowa forma eliminacji Gaussa . Komputery zwykle rozwiązują kwadratowe układy równań liniowych za pomocą rozkładu LU, a jest to również kluczowy krok podczas odwracania macierzy lub obliczania wyznacznika macierzy. Rozkład LU został wprowadzony przez polskiego matematyka Tadeusza Banachiewicza w 1938 roku.
Definicje
Niech A będzie macierzą kwadratową. LU na czynniki dotyczy faktoryzacji A , wraz z odpowiednim rzędzie i / lub uporządkowania kolumn lub permutacji, dla dwóch czynników - dolna matryca trójkątna L i górnej trójkątnej matrycy U :
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 3 × 3 A jej rozkład LU wygląda tak:
Bez odpowiedniego uporządkowania lub permutacji w macierzy faktoryzacja może się nie udać. Na przykład łatwo zweryfikować (poprzez rozwinięcie mnożenia macierzy), że . Jeśli , to co najmniej jeden z i musi być równy zero, co oznacza, że L lub U jest liczbą pojedynczą . Jest to niemożliwe, jeśli A jest nieosobliwe (odwracalne). To jest problem proceduralny. Można go usunąć po prostu zmieniając kolejność wierszy A tak, aby pierwszy element macierzy permutowanej był niezerowy. Ten sam problem w kolejnych krokach faktoryzacji można usunąć w ten sam sposób; zobacz podstawową procedurę poniżej.
Faktoryzacja LU z częściowym odchyleniem
Okazuje się, że odpowiednia permutacja wierszy (lub kolumn) jest wystarczająca do faktoryzacji LU. Faktoryzacja LU z częściowym odchyleniem (LUP) często odnosi się do faktoryzacji LU tylko z permutacjami wierszy:
gdzie L i U są znowu dolną i górną macierzą trójkątną, a P jest macierzą permutacji , która po pomnożeniu przez A , zmienia kolejność wierszy A . Okazuje się, że wszystkie macierze kwadratowe mogą być faktoryzowane 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 obrotem
Faktoryzacji LU pełnym wychyleniu obejmuje zarówno wierszy i kolumn permutacji:
gdzie L , U i P są zdefiniowane jak poprzednio, a Q jest macierzą permutacji, która zmienia kolejność kolumn A .
Rozkład dolna-przekątna-górna (LDU)
Dolna przekątnej, górne (LDU) rozkładu jest rozkład formie
gdzie D jest macierzą diagonalną , a L i U są macierzami unitargular , co oznacza, że wszystkie wpisy na przekątnych L i U są jednym.
Matryce prostokątne
Powyżej wymagaliśmy, aby A było macierzą kwadratową, ale 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 zawierają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 „wierszowa forma schodkowa” macierzy A .
Przykład
Faktoryzujemy następującą macierz 2 na 2:
Jednym ze sposobów znalezienia rozkładu LU tej prostej macierzy byłoby po prostu rozwiązanie równań liniowych przez inspekcję. Rozszerzenie mnożenia macierzy daje
Ten układ równań jest niedookreślony . W tym przypadku dowolne dwa niezerowe elementy macierzy L i U są parametrami rozwiązania i mogą być dowolnie ustawione 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 macierz trójkątna L była jednostkową macierzą trójkątną (tzn. ustawić wszystkie wpisy jej głównej przekątnej na jedynki). Wtedy układ równań ma następujące rozwiązanie:
Podstawiając te wartości do rozkładu LU powyżej uzysków
Istnienie i wyjątkowość
Macierze kwadratowe
Każda macierz kwadratowa dopuszcza faktoryzację LUP i PLU . Jeśli jest invertible , to dopuszcza faktoryzację LU (lub LDU ) wtedy i tylko wtedy, gdy wszystkie jego wiodące główne drugorzędne są niezerowe. Jeśli jest pojedynczą macierzą rangi , to dopuszcza faktoryzację LU, jeśli pierwsze wiodące główne drugorzędne są niezerowe, chociaż odwrotność nie jest prawdą.
Jeśli kwadratowa, odwracalna macierz ma LDU (faktoryzacja ze wszystkimi ukośnymi wpisami L i U równymi 1), to faktoryzacja jest unikalna. W takim przypadku faktoryzacja LU jest również unikalna, jeśli wymagamy, aby przekątna (lub ) składała się z jedynek.
Symetryczne macierze dodatnio określone
Jeśli jest symetryczna (lub hermitowskie jeśli jest kompleks) określoność formy , można rozmieścić tak, aby spraw U jest transpozycją sprzężoną z L . Oznacza to, że możemy napisać A jako
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. Co więcej, obliczanie rozkładu Cholesky'ego jest bardziej wydajne i numerycznie bardziej stabilne niż obliczanie niektórych innych rozkładów LU.
Macierze ogólne
W przypadku (niekoniecznie odwracalnej) macierzy nad dowolnym polem znane są dokładne warunki konieczne i wystarczające, w których ma ona faktoryzację LU. Warunki są wyrażone w postaci rang niektórych podmacierzy. Algorytm eliminacji Gaussa do uzyskiwania rozkładu LU również został rozszerzony do tego najbardziej ogólnego przypadku.
Algorytmy
Zamknięta formuła
Gdy faktoryzacja LDU istnieje i jest unikalna, istnieje zamknięta (wyraźna) formuła dla elementów L , D i U w kategoriach stosunków wyznaczników pewnych podmacierzy pierwotnej macierzy A . W szczególności , i dla , to stosunek -tej głównej podmacierzy do -tej głównej podmacierzy. Obliczenie wyznaczników jest obliczeniowo kosztowne , więc ta jednoznaczna formuła nie jest stosowana w praktyce.
Korzystanie z eliminacji Gaussa
Poniższy algorytm jest zasadniczo zmodyfikowaną formą eliminacji Gaussa . Obliczanie rozkładu LU przy użyciu tego algorytmu wymaga operacji zmiennoprzecinkowych, ignorując warunki niższego rzędu. Częściowe obracanie dodaje tylko wyraz kwadratowy; tak nie jest w przypadku pełnego obrotu.
Mając macierz N × N , zdefiniuj .
Eliminujemy elementy macierzy poniżej głównej przekątnej w n- tej kolumnie A ( n −1) dodając do i -tego wiersza tej macierzy n- ty wiersz pomnożony przez
Można to zrobić mnożąc A ( n − 1) w lewo przez dolną macierz trójkątną
czyli macierz jednostkowa N × N z n -tą kolumną zastąpioną wektorem
Ustawiamy
Po N − 1 krokach wyeliminowaliśmy wszystkie elementy macierzy poniżej głównej przekątnej, więc otrzymujemy górną macierz trójkątną A ( N − 1 ) . Znajdujemy rozkład
Oznacz górną macierz trójkątną A ( N − 1) przez U , i . Ze względu na odwrotną niższej trójkątnej matrycy L brak jest znowu trójkątny dolna matryca i mnożenie dwóch dolnych trójkątnych macierzy jest znowu mniejsza matryca trójkątna wynika, że L jest niższy matryca trójkątna. Co więcej, widać, że
Otrzymujemy
Oczywiste jest, że aby ten algorytm działał, trzeba mieć go na każdym kroku (patrz definicja ). Jeśli to założenie w którymś momencie się nie powiedzie, przed kontynuowaniem należy zamienić n -ty wiersz na inny wiersz poniżej. Dlatego generalnie wygląda dekompozycja LU .
Zauważ, że rozkład uzyskany w tej procedurze jest rozkładem Doolittle'a : główna przekątna L składa się wyłącznie z 1 sekundy. Jeśli przystąpimy do usuwania elementów powyżej głównej przekątnej dodając wielokrotności kolumn (zamiast usuwać elementy poniżej przekątnej dodając wielokrotności rzędów ), otrzymalibyśmy dekompozycję Crout , gdzie główna przekątna U wynosi 1 s .
Innym (równoważnym) sposobem wytwarzania dekompozycji Crouta danej macierzy A jest otrzymanie dekompozycji Doolittle'a transpozycji A . Rzeczywiście, jeśli jest to dekompozycja LU uzyskana za pomocą algorytmu przedstawionego w tej sekcji, to biorąc i , otrzymujemy dekompozycję Crouta.
Poprzez rekurencję
Cormen i in. opisać rekurencyjny algorytm dekompozycji LUP.
Mając macierz A , niech P 1 będzie macierzą permutacji taką, że
- ,
gdzie , jeśli istnieje niezerowy wpis w pierwszej kolumnie A ; lub podjąć P 1 jako macierz tożsamości inaczej. Teraz niech , jeśli ; lub inaczej. Mamy
Teraz możemy rekurencyjnie znaleźć rozkład LUP . Niech . W związku z tym
który jest rozkładem LUP A .
Algorytm randomizowany
Możliwe jest znalezienie aproksymacji niskiego rzędu do rozkładu LU przy użyciu algorytmu randomizowanego. Biorąc pod uwagę macierz wejściowego i pożądany niski stopień , gdy powraca LU permutacji randomizowanych matryc i górny / dolny trapezowe matryce o rozmiarze i odpowiednio, tak, że z dużym prawdopodobieństwem , jeżeli jest stałą zależną od parametrów algorytmu i jest -tym wartość osobliwa macierzy wejściowej .
Teoretyczna złożoność
Jeśli dwie macierze rzędu n można pomnożyć w czasie M ( n ), gdzie M ( n ) ≥ n a dla pewnego n > 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 .
Rozkład macierzowy rzadki
Opracowano specjalne algorytmy do faktoryzacji dużych macierzy rzadkich . Algorytmy te próbują znaleźć czynniki rzadkie L i U . Idealnie, koszt obliczeń jest określony 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 uporządkowań, które minimalizują wypełnienie, można rozwiązać za pomocą teorii grafów .
Aplikacje
Rozwiązywanie równań liniowych
Mając dany układ równań liniowych w postaci macierzowej
chcemy rozwiązać równanie dla x , biorąc pod uwagę A i b . Załóżmy, że uzyskaliśmy już rozkład LUP A taki, że , więc .
W tym przypadku rozwiązanie odbywa się w dwóch logicznych krokach:
- Najpierw rozwiązujemy równanie dla y .
- Po drugie, rozwiązujemy równanie na x .
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 i wstecz bez użycia procesu eliminacji Gaussa (jednak potrzebujemy tego procesu lub odpowiednika do obliczenia samego rozkładu LU ).
Powyższa procedura może być wielokrotnie stosowana do wielokrotnego rozwiązywania równania dla różnych b . W tym przypadku szybciej (i wygodniej) wykonać dekompozycję LU macierzy A raz, a następnie rozwiązać macierze trójkątne dla różnych b , zamiast za każdym razem stosować eliminację Gaussa. Można by sądzić, że macierze L i U „zakodowały” proces eliminacji Gaussa.
Koszt rozwiązania układu równań liniowych to w przybliżeniu operacje zmiennoprzecinkowe, jeśli macierz ma rozmiar . To sprawia, że jest dwa razy szybszy niż algorytmy oparte na dekompozycji QR , które kosztują operacje zmiennoprzecinkowe, gdy używane są odbicia Householder . Z tego powodu zazwyczaj 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 -by- p , więc próbujemy znaleźć macierz X (również macierz n -by- p ):
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ą jednostkową o rozmiarze n . Wynikałoby z tego, że wynik X musi być odwrotnością A .
Obliczanie wyznacznika
Biorąc pod uwagę rozkład LUP macierzy kwadratowej A , wyznacznik A można obliczyć bezpośrednio jako
Drugie równanie wynika z faktu, że wyznacznik macierzy trójkątnej jest po prostu iloczynem jej przekątnych, a wyznacznik macierzy permutacyjnej jest równy (−1) S gdzie S jest liczbą wymian wierszy w dekompozycji .
W przypadku dekompozycji LU z pełnym obrotem, równa się również prawej stronie powyższego równania, jeśli przyjmiemy, że S będzie całkowitą liczbą wymian wierszy i kolumn.
Ta sama metoda łatwo stosuje się do rozkładu LU, ustawiając P równe macierzy jednostkowej.
Przykłady kodu
Przykład kodu C
/* INPUT: A - array of pointers to rows of a square matrix having dimension N
* Tol - small tolerance number to detect failure when the matrix is near degenerate
* OUTPUT: Matrix A is changed, it contains a copy of both matrices L-E and U as A=(L-E)+U such that P*A=L*U.
* The permutation matrix is not stored as a matrix, but in an integer vector P of size N+1
* containing column indexes where the permutation matrix has "1". The last element P[N]=S+N,
* where S is the number of row exchanges needed for determinant computation, det(P)=(-1)^S
*/
int LUPDecompose(double **A, int N, double Tol, int *P) {
int i, j, k, imax;
double maxA, *ptr, absA;
for (i = 0; i <= N; i++)
P[i] = i; //Unit permutation matrix, P[N] initialized with N
for (i = 0; i < N; i++) {
maxA = 0.0;
imax = i;
for (k = i; k < N; k++)
if ((absA = fabs(A[k][i])) > maxA) {
maxA = absA;
imax = k;
}
if (maxA < Tol) return 0; //failure, matrix is degenerate
if (imax != i) {
//pivoting P
j = P[i];
P[i] = P[imax];
P[imax] = j;
//pivoting rows of A
ptr = A[i];
A[i] = A[imax];
A[imax] = ptr;
//counting pivots starting from N (for determinant)
P[N]++;
}
for (j = i + 1; j < N; j++) {
A[j][i] /= A[i][i];
for (k = i + 1; k < N; k++)
A[j][k] -= A[j][i] * A[i][k];
}
}
return 1; //decomposition done
}
/* INPUT: A,P filled in LUPDecompose; b - rhs vector; N - dimension
* OUTPUT: x - solution vector of A*x=b
*/
void LUPSolve(double **A, int *P, double *b, int N, double *x) {
for (int i = 0; i < N; i++) {
x[i] = b[P[i]];
for (int k = 0; k < i; k++)
x[i] -= A[i][k] * x[k];
}
for (int i = N - 1; i >= 0; i--) {
for (int k = i + 1; k < N; k++)
x[i] -= A[i][k] * x[k];
x[i] /= A[i][i];
}
}
/* INPUT: A,P filled in LUPDecompose; N - dimension
* OUTPUT: IA is the inverse of the initial matrix
*/
void LUPInvert(double **A, int *P, int N, double **IA) {
for (int j = 0; j < N; j++) {
for (int i = 0; i < N; i++) {
IA[i][j] = P[i] == j ? 1.0 : 0.0;
for (int k = 0; k < i; k++)
IA[i][j] -= A[i][k] * IA[k][j];
}
for (int i = N - 1; i >= 0; i--) {
for (int k = i + 1; k < N; k++)
IA[i][j] -= A[i][k] * IA[k][j];
IA[i][j] /= A[i][i];
}
}
}
/* INPUT: A,P filled in LUPDecompose; N - dimension.
* OUTPUT: Function returns the determinant of the initial matrix
*/
double LUPDeterminant(double **A, int *P, int N) {
double det = A[0][0];
for (int i = 1; i < N; i++)
det *= A[i][i];
return (P[N] - N) % 2 == 0 ? det : -det;
}
Przykład kodu C#
public class SystemOfLinearEquations
{
public double[] SolveUsingLU(double[,] matrix, double[] rightPart, int n)
{
// decomposition of matrix
double[,] lu = new double[n, n];
double sum = 0;
for (int i = 0; i < n; i++)
{
for (int j = i; j < n; j++)
{
sum = 0;
for (int k = 0; k < i; k++)
sum += lu[i, k] * lu[k, j];
lu[i, j] = matrix[i, j] - sum;
}
for (int j = i + 1; j < n; j++)
{
sum = 0;
for (int k = 0; k < i; k++)
sum += lu[j, k] * lu[k, i];
lu[j, i] = (1 / lu[i, i]) * (matrix[j, i] - sum);
}
}
// lu = L+U-I
// find solution of Ly = b
double[] y = new double[n];
for (int i = 0; i < n; i++)
{
sum = 0;
for (int k = 0; k < i; k++)
sum += lu[i, k] * y[k];
y[i] = rightPart[i] - sum;
}
// find solution of Ux = y
double[] x = new double[n];
for (int i = n - 1; i >= 0; i--)
{
sum = 0;
for (int k = i + 1; k < n; k++)
sum += lu[i, k] * x[k];
x[i] = (1 / lu[i, i]) * (y[i] - sum);
}
return x;
}
}
Przykład kodu MATLAB
function LU = LUDecompDoolittle(A)
n = length(A);
LU = A;
% decomposition of matrix, Doolittle’s Method
for i = 1:1:n
for j = 1:(i - 1)
LU(i,j) = (LU(i,j) - LU(i,1:(j - 1))*LU(1:(j - 1),j)) / LU(j,j);
end
j = i:n;
LU(i,j) = LU(i,j) - LU(i,1:(i - 1))*LU(1:(i - 1),j);
end
%LU = L+U-I
end
function x = SolveLinearSystem(LU, B)
n = length(LU);
y = zeros(size(B));
% find solution of Ly = B
for i = 1:n
y(i,:) = B(i,:) - L(i,1:i)*y(1:i,:);
end
% find solution of Ux = y
x = zeros(size(B));
for i = n:(-1):1
x(i,:) = (y(i,:) - U(i,(i + 1):n)*x((i + 1):n,:))/U(i, i);
end
end
A = [ 4 3 3; 6 3 3; 3 4 3 ]
LU = LUDecompDoolittle(A)
B = [ 1 2 3; 4 5 6; 7 8 9; 10 11 12 ]'
x = SolveLinearSystem(LU, B)
A * x
Zobacz też
- Rozkład bloku LU
- Rozkład Bruhata
- Rozkład Choleskiego
- Niepełna faktoryzacja LU
- Redukcja LU
- Rozkład macierzy
- Rozkład QR
Uwagi
Bibliografia
- Wiązka, James R.; Hopcroft, John (1974), „trójkątna faktoryzacja i odwracanie przez szybkie mnożenie macierzy”, Matematyka obliczeń , 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), Obliczenia macierzy (3rd ed.), 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 rozdział 3.5. N − 1
- Gospodarz, Alston S. (1975), Teoria macierzy w analizie numerycznej , New York: Dover Publications , MR 0378371.
- Okunev, Paweł; Johnson, Charles R. (1997), Warunki konieczne i wystarczające dla istnienia faktoryzacji LU macierzy arbitralnej , arXiv : math.NA/0506382.
- Poole, David (2006), Algebra Liniowa: Nowoczesne Wprowadzenie (2nd ed.), Kanada: Thomson Brooks / Cole, ISBN 978-0-534-99845-5.
- Prasa, WH; Teukolski SA; Vetterling, WT; Flannery, BP (2007), „Sekcja 2.3” , Przepisy numeryczne: The Art of Scientific Computing (3rd ed.), New York: 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.
Zewnętrzne linki
Bibliografia
- Rozkład LU w MathWorld .
- Rozkład LU na Math-Linuksie .
- Dekompozycja LU w Holistycznym Instytucie Metod Numerycznych
- Faktoryzacja macierzy LU . Odniesienie MATLAB.
Kod komputerowy
- LAPACK to zbiór podprogramów FORTRAN do rozwiązywania gęstych problemów algebry liniowej
- ALGLIB zawiera częściowy port LAPACK do C++, C#, Delphi itp.
- Kod C++ , prof. J. Loomis, University of Dayton
- Kod C , biblioteka źródeł matematyki
- LU w X10
Zasoby online
- WebApp opisowo rozwiązujący układy równań liniowych za pomocą dekompozycji LU
- Kalkulator macierzy , bluebit.gr
- Narzędzie do rozkładu LU , uni-bonn.de
- Rozkład LU autorstwa Eda Pegga, Jr. , Projekt demonstracji Wolfram , 2007.