Rozklad LU - LU decomposition

V numerické analýze a lineární algebře , nižší-horní ( LU ), rozklad nebo faktorizace faktorů A matrici jako produkt nižší trojúhelníkovou matici a horní trojúhelníkovou matici. Produkt někdy obsahuje také permutační matici . Rozklad LU lze považovat za maticovou formu Gaussovy eliminace . Počítače obvykle řeší čtvercové systémy lineárních rovnic pomocí dekompozice LU a je také klíčovým krokem při invertování matice nebo výpočtu determinantu matice. Rozklad LU zavedl polský matematik Tadeusz Banachiewicz v roce 1938.

Definice

LDU rozklad Walshovy matice

Nechť A je čtvercová matice. LU faktorizace odkazuje na faktorizace A , se správnou řádků a / nebo sloupců orderings nebo obměn, do dvou faktorů - nižší trojúhelníková matice L a horní trojúhelníkovou matici U :

${\ displaystyle A = LU.}$

V dolní trojúhelníkové matici jsou všechny prvky nad úhlopříčkou nulové, v horní trojúhelníkové matici jsou všechny prvky pod úhlopříčkou nulové. Například pro matici 3 × 3 A vypadá její rozklad LU takto:

${\ displaystyle {\ begin {bmatrix} a_ {11} & a_ {12} & a_ {13} \\ a_ {21} & a_ {22} & a_ {23} \\ a_ {31} & a_ {32} & 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 { bmatrix}}.}$

Bez řádného uspořádání nebo permutací v matici se faktorizace nemusí zhmotnit. Například je snadné ověřit (rozšířením maticového násobení), že . Pokud , pak alespoň jeden z a musí být nula, což znamená, že buď L nebo U je singulární . To je nemožné, pokud A je nesingulární (invertibilní). Jedná se o procesní problém. Lze jej odstranit jednoduchým přeuspořádáním řádků A tak, aby první prvek permutované matice byl nenulový. Stejný problém v následných faktorizačních krocích lze odstranit stejným způsobem; viz základní postup níže. ${\ textstyle a_ {11} = \ ell _ {11} u_ {11}}$${\ textstyle a_ {11} = 0}$${\ textstyle \ ell _ {11}}$${\ textstyle u_ {11}}$

Faktorizace LU s částečným otáčením

Ukazuje se, že pro faktorizaci LU stačí správná permutace v řádcích (nebo sloupcích). Faktorizace LU s částečným otáčením (LUP) často odkazuje na faktorizaci LU pouze s permutacemi řádků:

${\ displaystyle PA = LU,}$

kde L a U jsou opět spodní a horní trojúhelníkové matice, a P je permutace matice , která při ponechání násobený na A , přeřadí řádky A . Ukazuje se, že všechny čtvercové matice lze v této formě faktorizovat a faktorizace je v praxi numericky stabilní. Díky tomu je rozklad LUP užitečnou technikou v praxi.

Faktorizace LU s plným otáčením

LU faktorizace s plnou otáčení zahrnuje jak řádků a sloupců obměny:

${\ displaystyle PAQ = LU,}$

kde L , U a P mají výše uvedený význam, a Q je permutací matice, která mění pořadí sloupců A .

Rozklad dolní diagonální-horní (LDU)

Nižší diagonální-horní (LDU) rozklad je rozklad podobě

${\ displaystyle A = LDU,}$

kde D je diagonální matice a L a U jsou unitriangulární matice , což znamená, že všechny položky na diagonálech L a U jsou jedna.

Obdélníkové matice

Výše jsme požadovali, aby A byla čtvercová matice, ale tyto dekompozice lze také zobecnit na obdélníkové matice. V tom případě, L a D jsou čtvercové matice, které oba mají stejný počet řádků jako A , a U má přesně stejné rozměry jako A . Horní trojúhelníkový by měl být interpretován tak, že má pouze nulové záznamy pod hlavní úhlopříčkou, která začíná v levém horním rohu. Stejně tak přesnější termín pro U je, že se jedná o „řádek sledu forma“ matice A .

Příklad

Faktorizujeme následující matici 2 na 2:

${\ displaystyle {\ begin {bmatrix} 4 & 3 \\ 6 & 3 \ end {bmatrix}} = {\ begin {bmatrix} \ ell _ {11} & 0 \\\ ell _ {21} & \ ell _ {22} \ end {bmatrix}} {\ begin {bmatrix} u_ {11} & u_ {12} \\ 0 & u_ {22} \ end {bmatrix}}.}$

Jedním ze způsobů, jak najít rozklad LU této jednoduché matice, by bylo jednoduše vyřešit lineární rovnice kontrolou. Rozšiřování násobení matice dává

${\ Displaystyle {\ begin {aligned} \ 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 {aligned}}}$

Tento systém rovnic je poddeterminovaný . V tomto případě jsou libovolné dva nenulové prvky matic L a U parametry řešení a lze je libovolně nastavit na libovolnou nenulovou hodnotu. Abychom našli jedinečný rozklad LU, je nutné dát určité omezení na matice L a U. Například můžeme pohodlně požadovat, aby spodní trojúhelníková matice L byla jednotková trojúhelníková matice (tj. Nastavila všechny položky její hlavní úhlopříčky na jedničky). Pak má soustava rovnic následující řešení:

${\ Displaystyle {\ begin {aligned} \ ell _ {11} = \ ell _ {22} & = 1 \\\ ell _ {21} & = 1,5 \\ u_ {11} & = 4 \\ u_ {12 } & = 3 \\ u_ {22} & =-1,5 \ konec {zarovnáno}}}$

Substituce těchto hodnot do LU rozkladu nad výtěžky

${\ Displaystyle {\ begin {bmatrix} 4 & 3 \\ 6 & 3 \ end {bmatrix}} = {\ begin {bmatrix} 1 & 0 \\ 1,5 & 1 \ end {bmatrix}} {{begin {bmatrix} 4 & 3 \\ 0 & -1,5 \ end {bmatrix}}.}$

Existence a jedinečnost

Čtvercové matice

Jakákoli čtvercová matice připouští faktorizace LUP a PLU . Pokud je nevratný , pak připouští faktorizaci LU (nebo LDU ) právě tehdy, pokud jsou všichni jeho hlavní hlavní nezletilí nenulové. Pokud je singulární matice hodnosti , pak připouští faktorizaci LU, pokud jsou první přední hlavní nezletilí nenulové, ačkoli opak není pravdivý. ${\ textstyle A}$${\ textstyle A}$${\ textstyle A}$${\ textstyle k}$${\ textstyle k}$

Pokud má čtvercová, invertibilní matice LDU (faktorizace se všemi diagonálními vstupy L a U rovnou 1), pak je faktorizace jedinečná. V takovém případě je faktorizace LU také jedinečná, pokud požadujeme, aby úhlopříčka (nebo ) sestávala z jedniček. ${\ textstyle L}$${\ textstyle U}$

Symetrické matice s jednoznačnou pozitivitou

Pokud je symetrický (nebo Hermitian , pokud je komplex) pozitivně definitní matice , můžeme zajistit věci, takže U je konjugovat přemístit z L . To znamená, že můžeme napsat A jako

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

Tento rozklad se nazývá Choleskyho rozklad . Rozklad Cholesky vždy existuje a je jedinečný - za předpokladu, že je matice kladně definitivní. Kromě toho je výpočet Choleskyho dekompozice efektivnější a numericky stabilnější než výpočet některých jiných dekompozic LU.

Obecné matice

Pro (ne nutně invertibilní) matici v jakémkoli poli jsou známy přesné nezbytné a dostatečné podmínky, za kterých má LU faktorizaci. Podmínky jsou vyjádřeny v řadách určitých dílčích matic. Na tento nejobecnější případ byl také rozšířen Gaussův eliminační algoritmus pro získání dekompozice LU.

Algoritmy

Uzavřený vzorec

Pokud existuje LDU faktorizace a je jedinečný, je uzavřený (explicitní) vzorec pro prvky L , D , a U , pokud jde o poměry determinant některých submatic původní matice A . Zejména, and for , je poměr -th hlavní submatice k -th hlavní submatice. Výpočet determinantů je výpočetně nákladný , proto se tento explicitní vzorec v praxi nepoužívá. ${\ textstyle D_ {1} = A_ {1,1}}$${\ textstyle i = 2, \ ldots, n}$${\ textstyle D_ {i}}$${\ textstyle i}$${\ textstyle (i-1)}$

Pomocí Gaussovy eliminace

Následující algoritmus je v podstatě upravenou formou Gaussovy eliminace . Výpočet dekompozice LU pomocí tohoto algoritmu vyžaduje operace s plovoucí desetinnou čárkou, ignorování výrazů nižšího řádu. Částečné otáčení přidává pouze kvadratický výraz; toto neplatí pro úplné otočení. ${\ displaystyle {\ tfrac {2} {3}} n^{3}}$

Vzhledem k matici N × N definujte . ${\ Displaystyle A = (a_ {i, j}) _ {1 \ leq i, j \ leq N}}$${\ displaystyle A^{(0)}: = A}$

Eliminujeme maticové prvky pod hlavní úhlopříčkou v n -tém sloupci A ^{( n −1)} tak, že do i -tého řádku této matice přičteme n -tý řádek vynásobený

${\ Displaystyle-\ ell _ {i, n}: =-{\ frac {a_ {i, n}^{(n-1)}} {a_ {n, n}^{(n-1)}} }, \ quad i = n+1, \ dotsc, N.}$

To lze provést vynásobením A ^{( n  - 1)} vlevo spodní trojúhelníkovou maticí

${\ displaystyle L_ {n} = {\ begin {pmatrix} 1 &&&&& \ \ & \ ddots &&&&& \\ && 1 &&& \ \ & &&-\ ell _ {n+1, n} &&& \\ && \ vdots && \ ddots & \\ &&-\ ell _ {N, n} &&& 1 \ end {pmatrix}},}$

tj. matice identity N × N s jejím n -tým sloupcem nahrazeným vektorem${\ Displaystyle {\ begin {pmatrix} 0 & \ dotsm & 0 & 1 &-\ ell _ {n+1, n} & \ dotsm &-\ ell _ {N, n} \ end {pmatrix}}^{\ textyf {T} }.}$

Jsme si stanovili

${\ Displaystyle A^{(n)}: = L_ {n} A^{(n-1)}.}$

Po N  - 1 krocích jsme eliminovali všechny maticové prvky pod hlavní úhlopříčkou, takže získáme horní trojúhelníkovou matici A ^{( N  - 1)} . Nacházíme rozklad

${\ 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)}.}$

Horní trojúhelníkovou matici A ^{( N  - 1)} označte U , a . Protože inverze nižší trojúhelníkové matice L _n je opět nižší trojúhelníková matice a násobení dvou nižších trojúhelníkových matic je opět nižší trojúhelníková matice, vyplývá z toho, že L je nižší trojúhelníková matice. Navíc je to vidět ${\ textstyle L = L_ {1}^{-1} \ dotsm L_ {N-1}^{-1}}$

${\ displaystyle L = {\ begin {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 {pmatrix}}.}$

Získáváme ${\ displaystyle A = LU.}$

Je jasné, že aby tento algoritmus fungoval, musí mít každý krok (viz definice ). Pokud tento předpoklad v určitém okamžiku selže, je třeba před pokračováním vyměnit n -tý řádek s dalším řádkem pod ním. Z tohoto důvodu vypadá rozklad LU obecně . ${\ Displaystyle a_ {n, n}^{(n-1)} \ neq 0}$${\ displaystyle \ ell _ {i, n}}$${\ displaystyle P^{-1} A = LU}$

Rozklad získaný tímto postupem je Doolittlův rozklad : hlavní úhlopříčka L se skládá pouze z 1 s. Pokud bychom postupovali odstraněním prvků nad hlavní úhlopříčkou přidáním násobků sloupců (namísto odebírání prvků pod úhlopříčkou přidáním násobků řádků ), získali bychom Croutův rozklad , kde hlavní úhlopříčka U je 1 s .

Další (ekvivalent) způsob výroby Crout rozklad dané matice A je získat Doolittle rozklad přemístit A . Pokud je dekompozice LU získána pomocí algoritmu uvedeného v této části, pak pomocí a máme, že je to Croutův rozklad. ${\ textstyle A^{\ textyf {T}} = L_ {0} U_ {0}}$${\ textstyle L = U_ {0}^{\ textyf {T}}}$${\ textstyle U = L_ {0}^{\ textyf {T}}}$${\ displaystyle A = LU}$

Prostřednictvím rekurze

Cormen a kol. popsat rekurzivní algoritmus pro rozklad LUP.

Vzhledem k matici A nechť P ₁ je taková permutační matice

${\ Displaystyle P_ {1} A = \ left ({\ begin {array} {c | ccc} a && w^{\ textf {T}} & \\\ hline &&&& \\ v && A '& \\ &&& \ end {array }}\že jo)}$,

kde , pokud je v prvním sloupci A nenulový záznam ; nebo jinak vezměte P ₁ jako matici identity. Nyní nechť , pokud ; nebo jinak. My máme ${\ textstyle a \ neq 0}$${\ textstyle c = 1/a}$${\ textstyle a \ neq 0}$${\ textstyle c = 0}$

${\ displaystyle P_ {1} A = \ left ({\ begin {array} {c | ccc} 1 && 0 & \\\ hline &&& \ \ cv && I_ {n-1} & \\ &&& \ end {array}} \ right) \ left ({\ begin {array} {c | c} a & w^{\ textyf {T}} \\\ hline & \\ 0 & A'-cvw^{\ textyf {T}} \\ & \ end {pole} }\že jo).}$

Nyní můžeme rekurzivně najít rozklad LUP . Nech . Proto ${\ textstyle P '\ left (A'-cvw^{\ textf {T}} \ right) = L'U'}$${\ textstyle v '= P'v}$

${\ displaystyle \ left ({\ begin {array} {c | ccc} 1 && 0 & \\\ hline &&&& \\ 0 && P '& \\ &&& \ end {array}} \ right) P_ {1} A = \ left ({ \ begin {array} {c | ccc} 1 && 0 & \\\ hline &&& \ \ cv '&& L' & \\ &&& \ end {array}} \ right) \ left ({\ begin {array} {c | ccc} a && w ^{\ textyf {T}} & \\\ hline &&& \\ 0 && U '& \\ &&&& \ end {array}} \ right),}$

což je LUP rozklad A .

Randomizovaný algoritmus

Je možné najít aproximaci nízké hodnosti k rozkladu LU pomocí randomizovaného algoritmu. Vzhledem k tomu, vstupní matrici a požadovanou nízkou hodnost , randomizované LU vrátí permutační matrice a spodní / horní trapézové matice o velikosti a příslušně, tak, že s velkou pravděpodobností , kde je konstanta, která závisí na parametrech algoritmu a je tého singulární hodnota vstupní matice . ${\ textstyle A}$${\ textstyle k}$${\ textstyle P, Q}$${\ textstyle L, U}$${\ textstyle m \ times k}$${\ textstyle k \ times n}$${\ textstyle \ left \ | PAQ-LU \ right \ | _ {2} \ leq C \ sigma _ {k+1}}$${\ textstyle C}$${\ textstyle \ sigma _ {k+1}}$${\ textstyle (k+1)}$${\ textstyle A}$

Teoretická náročnost

Pokud lze v čase M ( n ) vynásobit dvě matice řádu n , kde M ( n ) ≥ n ^a pro nějaké n > 2, pak lze rozklad LU vypočítat v čase O ( M ( n )). To například znamená, že existuje algoritmus O ( n ^2,376 ) založený na algoritmu Coppersmith – Winograd .

Rozklad řídké matice

Byly vyvinuty speciální algoritmy pro faktorizaci velkých řídkých matic . Tyto algoritmy pokoušejí najít řídké faktory L a U . V ideálním případě jsou náklady na výpočet určeny spíše počtem nenulových záznamů než velikostí matice.

Tyto algoritmy využívají svobodu pro výměnu řádků a sloupců k minimalizaci vyplňování (položky, které se během provádění algoritmu mění z počáteční nuly na nenulovou hodnotu).

Obecné zacházení s objednávkami, které minimalizují vyplňování, lze řešit pomocí teorie grafů .

Aplikace

Řešení lineárních rovnic

Daný soustavou lineárních rovnic ve formě matice

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

chceme vyřešit rovnici pro x , danou A a b . Předpokládejme, že jsme již získali LUP rozklad A tak, že , takže . ${\ textstyle PA = LU}$${\ textstyle LU \ mathbf {x} = P \ mathbf {b}}$

V tomto případě se řešení provádí ve dvou logických krocích:

1. Nejprve vyřešíme rovnici pro y .${\ textstyle L \ mathbf {y} = P \ mathbf {b}}$
2. Za druhé, řešíme rovnici pro x .${\ textstyle U \ mathbf {x} = \ mathbf {y}}$

V obou případech máme co do činění s trojúhelníkovými maticemi ( L a U ), které lze vyřešit přímo dopřednou a zpětnou substitucí bez použití Gaussova eliminačního procesu (nicméně tento proces nebo ekvivalent potřebujeme k výpočtu samotného rozkladu LU ).

Výše uvedený postup lze opakovaně použít k řešení rovnice vícekrát pro různá b . V tomto případě je rychlejší (a pohodlnější) provést LU rozklad matice A jednou a poté vyřešit trojúhelníkové matice pro různá b , namísto použití Gaussovy eliminace pokaždé. O matricích L a U by se dalo předpokládat, že „zakódovaly“ Gaussův eliminační proces.

Náklady na řešení soustavy lineárních rovnic jsou přibližně operace s pohyblivou řádovou čárkou, pokud má matice velikost . Díky tomu je dvakrát rychlejší než algoritmy založené na dekompozici QR , která stojí za použití operací s plovoucí desetinnou čárkou při použití odrazů Householder . Z tohoto důvodu je obvykle preferován rozklad LU. ${\ textstyle {\ frac {2} {3}} n^{3}}$${\ textstyle A}$${\ textstyle n}$${\ textstyle {\ frac {4} {3}} n^{3}}$

Invertování matice

Při řešení soustavy rovnic, b je obvykle považován za vektor s délkou, která se rovná výšce matice A . V inverzi matice však místo vektoru b máme matici B , kde B je matice n -by -p , takže se snažíme najít matici X (také matici n -by -p ):

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

Pro řešení pro každý sloupec matice X můžeme použít stejný algoritmus, který byl uveden dříve . Nyní předpokládejme, že B je matice identity velikosti n . Vyplývalo by z toho, že výsledek X musí být inverzní A .

Výpočet determinantu

Vzhledem k rozkladu LUP čtvercové matice A lze determinant A vypočítat přímo jako ${\ textstyle A = P^{-1} LU}$

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

Druhá rovnice vyplývá ze skutečnosti, že determinant trojúhelníkové matice je jednoduše součin jejích diagonálních záznamů a že determinant permutační matice je roven (−1) ^S, kde S je počet výměn řádků při rozkladu .

V případě dekompozice LU s plným otáčením se také rovná pravé straně výše uvedené rovnice, pokud necháme S být celkový počet výměn řádků a sloupců. ${\ textstyle \ det (A)}$

Stejná metoda snadno platí pro rozklad LU nastavením P rovným matici identity.

Příklady kódu

Příklad kódu 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;
}

Příklad kódu 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;
    }
}

Příklad kódu 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

Viz také

• Blokovat rozklad LU
• Bruhatův rozklad
• Choleskyho rozklad
• Neúplná faktorizace LU
• Redukce LU
• Rozklad matice
• Rozklad QR

Poznámky

Reference

• Parta, James R .; Hopcroft, John (1974), „Triangular factorization and inversion by fast matrix multiplication“, 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), Úvod do algoritmů , MIT Press a McGraw-Hill, ISBN 978-0-262-03293-3.
• Golub, Gene H .; Van Loan, Charles F. (1996), Matrix Computations (3. vyd.), Baltimore: Johns Hopkins, ISBN 978-0-8018-5414-9.
• Horn, Roger A .; Johnson, Charles R. (1985), Matrix Analysis , Cambridge University Press, ISBN 978-0-521-38632-6. Viz část 3.5. N  - 1
• Householder, Alston S. (1975), Theory of Matrices in Numerical Analysis , New York: Dover Publications , MR  0378371.
• Okunev, Pavel; Johnson, Charles R. (1997), Nezbytné a dostatečné podmínky pro existenci LU faktorizace libovolné matice , arXiv : math.NA/0506382.
• Poole, David (2006), Linear Algebra: A Modern Introduction (2. vyd.), Kanada: Thomson Brooks/Cole, ISBN 978-0-534-99845-5.
• Stiskněte, WH; Teukolsky, SA; Vetterling, WT; Flannery, BP (2007), „Oddíl 2.3“ , Numerical Recipes: The Art of Scientific Computing (3. vyd.), New York: Cambridge University Press, ISBN 978-0-521-88068-8.
• Trefethen, Lloyd N .; Bau, David (1997), Numerická lineární algebra , Philadelphia: Společnost pro průmyslovou a aplikovanou matematiku, ISBN 978-0-89871-361-9.

externí odkazy

Reference

• Rozklad LU na MathWorld .
• Rozklad LU na Math-Linux .
• Rozklad LU v Holistic Numerical Methods Institute
• Faktorizace matice LU . Reference MATLAB.

Počítačový kód

• LAPACK je sbírka podprogramů FORTRAN pro řešení hustých problémů lineární algebry
• ALGLIB obsahuje částečný port LAPACK do C ++, C#, Delphi atd.
• C ++ kód , Prof. J. Loomis, University of Dayton
• C kód , knihovna zdrojů matematiky
• LU v X10

Online zdroje

• WebApp popisně řeší systémy lineárních rovnic pomocí LU Decomposition
• Matrix Calculator , bluebit.gr
• Nástroj pro dekompozici LU , uni-bonn.de
• LU Decomposition od Ed Pegg, Jr. , The Wolfram Demonstrations Project , 2007.