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

Rozkład LDU macierzy Walsha

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 Umacierzami 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:

  1. Najpierw rozwiązujemy równanie dla y .
  2. 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ż

Uwagi

Bibliografia

Zewnętrzne linki

Bibliografia

Kod komputerowy

Zasoby online