Decomposição LU - LU decomposition

Na análise numérica e na álgebra linear , a decomposição inferior-superior ( LU ) ou fatoração fatores uma matriz como o produto de uma matriz triangular inferior e uma matriz triangular superior. O produto às vezes também inclui uma matriz de permutação . A decomposição LU pode ser vista como a forma de matriz de eliminação gaussiana . Os computadores geralmente resolvem sistemas quadrados de equações lineares usando a decomposição LU, e também é uma etapa fundamental ao inverter uma matriz ou calcular o determinante de uma matriz. A decomposição LU foi introduzida pelo matemático polonês Tadeusz Banachiewicz em 1938.

Definições

Decomposição LDU de uma matriz Walsh

Seja A uma matriz quadrada. Uma fatoração LU refere-se à fatoração de A , com ordenações ou permutações adequadas de linhas e / ou colunas, em dois fatores - uma matriz triangular inferior L e uma matriz triangular superior U :

${\ displaystyle A = LU.}$

Na matriz triangular inferior, todos os elementos acima da diagonal são zero; na matriz triangular superior, todos os elementos abaixo da diagonal são zero. Por exemplo, para uma matriz A 3 × 3 , sua decomposição LU é assim:

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

Sem uma ordenação adequada ou permutações na matriz, a fatoração pode não se materializar. Por exemplo, é fácil verificar (expandindo a multiplicação da matriz) isso . Se , então, pelo menos um de e tem que ser zero, o que implica que L ou U são singulares . Isso é impossível se A for não singular (invertível). Este é um problema processual. Ele pode ser removido simplesmente reordenando as linhas de A de modo que o primeiro elemento da matriz permutada seja diferente de zero. O mesmo problema nas etapas de fatoração subsequentes pode ser removido da mesma maneira; veja o procedimento básico abaixo. ${\ textstyle a_ {11} = \ ell _ {11} u_ {11}}$${\ textstyle a_ {11} = 0}$${\ textstyle \ ell _ {11}}$${\ textstyle u_ {11}}$

Fatoração LU com pivotamento parcial

Acontece que uma permutação adequada em linhas (ou colunas) é suficiente para a fatoração LU. Fatoração LU com pivotamento parcial (LUP) se refere frequentemente à fatoração LU apenas com permutações de linha:

${\ displaystyle PA = LU,}$

onde L e L são novamente inferior e superior matrizes triangulares, e P é uma matriz de permutação , que, quando deixou-multiplicado a um , reordena as linhas de Uma . Acontece que todas as matrizes quadradas podem ser fatoradas dessa forma, e a fatoração é numericamente estável na prática. Isso torna a decomposição do LUP uma técnica útil na prática.

Fatoração LU com pivotamento completo

Uma fatoração LU com pivotamento completo envolve permutações de linha e coluna:

${\ displaystyle PAQ = LU,}$

onde L , L e P são definidos como anteriormente, e Q é uma matriz de permutação que reordena as colunas de um .

Decomposição diagonal inferior-superior (LDU)

Uma decomposição inferior diagonal superior (LDU) é uma decomposição da forma

${\ displaystyle A = LDU,}$

onde D é uma matriz diagonal e L e U são matrizes unitriangulares , o que significa que todas as entradas nas diagonais de L e U são uma.

Matrizes retangulares

Acima exigimos que A seja uma matriz quadrada, mas essas decomposições podem ser generalizadas para matrizes retangulares também. Nesse caso, G e D são matrizes quadradas sendo que ambos têm o mesmo número de filas como um , e L possui exactamente as mesmas dimensões que um . O triangular superior deve ser interpretado como tendo apenas zero entradas abaixo da diagonal principal, que começa no canto superior esquerdo. Da mesma forma, o termo mais preciso para U é que é a "forma escalonada" da matriz A .

Exemplo

Fatoramos a seguinte matriz 2 por 2:

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

Uma maneira de encontrar a decomposição LU dessa matriz simples seria simplesmente resolver as equações lineares por inspeção. Expandir a multiplicação da matriz dá

${\ displaystyle {\ begin {alinhado} \ 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 {alinhado}}}$

Este sistema de equações é subdeterminado . Nesse caso, quaisquer dois elementos diferentes de zero das matrizes L e U são parâmetros da solução e podem ser definidos arbitrariamente para qualquer valor diferente de zero. Portanto, para encontrar a decomposição LU única, é necessário colocar algumas restrições nas matrizes L e U. Por exemplo, podemos convenientemente exigir que a matriz triangular inferior L seja uma matriz triangular unitária (ou seja, defina todas as entradas de sua diagonal principal como uns). Então, o sistema de equações tem a seguinte solução:

${\ displaystyle {\ begin {alinhado} \ ell _ {11} = \ ell _ {22} & = 1 \\\ ell _ {21} & = 1,5 \\ u_ {11} & = 4 \\ u_ {12 } & = 3 \\ u_ {22} & = - 1,5 \ fim {alinhado}}}$

Substituir esses valores na decomposição LU acima produz

${\ displaystyle {\ begin {bmatrix} 4 e 3 \\ 6 & 3 \ end {bmatrix}} = {\ begin {bmatrix} 1 & 0 \\ 1.5 & 1 \ end {bmatrix}} {\ begin {bmatrix} 4 e 3 \\ 0 & -1.5 \ fim de {bmatrix}}.}$

Existência e singularidade

Matrizes quadradas

Qualquer matriz quadrada admite fatorações LUP e PLU . Se for invertível , então admite uma fatoração LU (ou LDU ) se e somente se todos os seus principais menores principais forem diferentes de zero. Se for uma matriz de classificação singular , então ela admite uma fatoração LU se os primeiros principais menores forem diferentes de zero, embora o inverso não seja verdadeiro. ${\ textstyle A}$${\ textstyle A}$${\ textstyle A}$${\ textstyle k}$${\ textstyle k}$

Se uma matriz quadrada e invertível tem uma LDU (fatoração com todas as entradas diagonais de L e U iguais a 1), então a fatoração é única. Nesse caso, a fatoração LU também é única se exigirmos que a diagonal de (ou ) consista em uns. ${\ textstyle L}$${\ textstyle U}$

Matrizes simétricas positivas-definidas

Se A é uma simétrica (ou Hermitian , se A é complexa) matriz definida positiva , podemos organizar as coisas de modo que U é a transposta conjugada de L . Ou seja, podemos escrever A como

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

Essa decomposição é chamada de decomposição de Cholesky . A decomposição de Cholesky sempre existe e é única - desde que a matriz seja definida positiva. Além disso, o cálculo da decomposição de Cholesky é mais eficiente e numericamente mais estável do que o cálculo de algumas outras decomposições LU.

Matrizes gerais

Para uma matriz (não necessariamente invertível) em qualquer campo, as condições exatas necessárias e suficientes sob as quais ela tem uma fatoração LU são conhecidas. As condições são expressas em termos de classificação de certas submatrizes. O algoritmo de eliminação gaussiana para obter a decomposição LU também foi estendido para este caso mais geral.

Algoritmos

Fórmula fechada

Quando existe um fatoração LDU e é única, não existe uma fórmula fechada (explícita) para os elementos de L , D e L em termos de razões de determinantes de certos submatrizes da matriz original Uma . Em particular ,, e para , é a razão entre a -ésima submatriz principal e a -ésima submatriz principal. O cálculo dos determinantes é caro do ponto de vista computacional , portanto, essa fórmula explícita não é usada na prática. ${\ textstyle D_ {1} = A_ {1,1}}$${\ textstyle i = 2, \ ldots, n}$${\ textstyle D_ {i}}$${\ textstyle i}$${\ textstyle (i-1)}$

Usando a eliminação de Gauss

O algoritmo a seguir é essencialmente uma forma modificada de eliminação gaussiana . O cálculo de uma decomposição de LU usando este algoritmo requer operações de ponto flutuante, ignorando os termos de ordem inferior. O pivotamento parcial adiciona apenas um termo quadrático; este não é o caso do pivotamento completo. ${\ displaystyle {\ tfrac {2} {3}} n ^ {3}}$

Dado um N × N matriz , definir . ${\ displaystyle A = (a_ {i, j}) _ {1 \ leq i, j \ leq N}}$${\ displaystyle A ^ {(0)}: = A}$

Eliminamos os elementos da matriz abaixo da diagonal principal na n- ésima coluna de A ^{( n −1)} adicionando à i- ésima linha desta matriz a n- ésima linha multiplicada por

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

Isso pode ser feito multiplicando A ^{( n  - 1)} à esquerda com a matriz triangular inferior

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

isto é, o N x N matriz de identidade com a sua n coluna -ésimo substituída pelo vector${\ displaystyle {\ begin {pmatrix} 0 & \ dotsm & 0 & 1 & - \ ell _ {n + 1, n} & \ dotsm & - \ ell _ {N, n} \ end {pmatrix}} ^ {\ textf {T} }.}$

Montamos

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

Após N  - 1 passos, eliminamos todos os elementos da matriz abaixo da diagonal principal, de modo que obtemos uma matriz triangular superior A ^{( N  - 1)} . Encontramos a decomposição

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

Denote a matriz triangular superior A ^{( N  - 1)} por U , e . Como o inverso de uma matriz triangular inferior L _n é novamente uma matriz triangular inferior, e a multiplicação de duas matrizes triangulares inferiores é novamente uma matriz triangular inferior, segue-se que L é uma matriz triangular inferior. Além disso, pode-se ver que ${\ 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}}.}$

Nós obtemos ${\ displaystyle A = LU.}$

É claro que para que este algoritmo funcione, é necessário ter em cada etapa (ver a definição de ). Se essa suposição falhar em algum ponto, será necessário intercambiar a n- ésima linha com outra linha abaixo dela antes de continuar. É por isso que uma decomposição LU em geral se parece . ${\ displaystyle a_ {n, n} ^ {(n-1)} \ neq 0}$${\ displaystyle \ ell _ {i, n}}$${\ displaystyle P ^ {- 1} A = LU}$

Observe que a decomposição obtida por meio desse procedimento é uma decomposição Doolittle : a diagonal principal de L é composta apenas por 1 s. Se prosseguirmos removendo elementos acima da diagonal principal adicionando múltiplos das colunas (em vez de remover elementos abaixo da diagonal adicionando múltiplos das linhas ), obteríamos uma decomposição de Crout , onde a diagonal principal de U é de 1 s .

Outra forma (equivalente) de produzir uma decomposição Crout de uma dada matriz Uma é o de obter uma decomposição Doolittle da transposição de uma . De fato, se é a decomposição LU obtida através do algoritmo apresentado nesta seção, então tomando e , temos que é uma decomposição de Crout. ${\ textstyle A ^ {\ textf {T}} = L_ {0} U_ {0}}$${\ textstyle L = U_ {0} ^ {\ textf {T}}}$${\ textstyle U = L_ {0} ^ {\ textf {T}}}$${\ displaystyle A = LU}$

Por recursão

Cormen et al. descreve um algoritmo recursivo para decomposição de LUP.

Dada uma matriz A , seja P ₁ uma matriz de permutação tal que

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

onde , se houver uma entrada diferente de zero na primeira coluna de A ; ou tome P ₁ como a matriz de identidade de outra forma. Agora vamos , se ; ou de outra forma. Nós temos ${\ 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 ^ {\ textf {T}} \\\ hline & \\ 0 & A'-cvw ^ {\ textf {T}} \\ & \ end {array} }\direito).}$

Agora podemos encontrar recursivamente uma decomposição LUP . Deixe . Portanto ${\ 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 ^ {\ textf {T}} & \\\ hline &&& \\ 0 && U '& \\ &&& \ end {array}} \ right),}$

que é uma decomposição LUP de um .

Algoritmo Randomizado

É possível encontrar uma aproximação de classificação baixa para uma decomposição LU usando um algoritmo aleatório. Dada uma matriz de entrada e uma baixa classificação desejada , a LU aleatória retorna matrizes de permutação e matrizes trapezoidais inferior / superior de tamanho e , respectivamente, de modo que com alta probabilidade , onde é uma constante que depende dos parâmetros do algoritmo e é o -ésimo valor singular da matriz de entrada . ${\ 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}$

Complexidade teórica

Se duas matrizes de ordem n podem ser multiplicadas no tempo M ( n ), onde M ( n ) ≥ n ^a para algum n > 2, então uma decomposição LU pode ser calculada no tempo O ( M ( n )). Isso significa, por exemplo, que existe um algoritmo O ( n ^2.376 ) baseado no algoritmo Coppersmith – Winograd .

Decomposição de matriz esparsa

Algoritmos especiais foram desenvolvidos para fatorar grandes matrizes esparsas . Esses algoritmos tentam encontrar fatores esparsos L e U . Idealmente, o custo de computação é determinado pelo número de entradas diferentes de zero, e não pelo tamanho da matriz.

Esses algoritmos usam a liberdade de trocar linhas e colunas para minimizar o preenchimento (entradas que mudam de um zero inicial para um valor diferente de zero durante a execução de um algoritmo).

O tratamento geral de pedidos que minimizam o preenchimento pode ser feito usando a teoria dos grafos .

Formulários

Resolvendo equações lineares

Dado um sistema de equações lineares em forma de matriz

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

queremos resolver a equação para x , dados A e b . Suponha que já tenhamos obtido a decomposição LUP de A tal que , então . ${\ textstyle PA = LU}$${\ textstyle LU \ mathbf {x} = P \ mathbf {b}}$

Neste caso, a solução é feita em duas etapas lógicas:

1. Primeiro, resolvemos a equação para y .${\ textstyle L \ mathbf {y} = P \ mathbf {b}}$
2. Em segundo lugar, resolvemos a equação para x .${\ textstyle U \ mathbf {x} = \ mathbf {y}}$

Em ambos os casos, estamos lidando com matrizes triangulares ( L e U ), que podem ser resolvidas diretamente por substituição direta e regressiva sem usar o processo de eliminação gaussiana (no entanto, precisamos desse processo ou equivalente para calcular a decomposição LU em si).

O procedimento acima pode ser aplicado repetidamente para resolver a equação várias vezes para diferentes b . Nesse caso, é mais rápido (e mais conveniente) fazer uma decomposição LU da matriz A uma vez e, em seguida, resolver as matrizes triangulares para os diferentes b , em vez de usar a eliminação de Gauss a cada vez. Pode-se pensar que as matrizes L e U "codificaram" o processo de eliminação gaussiana.

O custo de resolver um sistema de equações lineares é aproximadamente operações de ponto flutuante se a matriz tiver tamanho . Isso o torna duas vezes mais rápido que algoritmos baseados em decomposição QR , que custa cerca de operações de ponto flutuante quando as reflexões de Householder são usadas. Por esta razão, a decomposição LU é geralmente preferida. ${\ textstyle {\ frac {2} {3}} n ^ {3}}$${\ textstyle A}$${\ textstyle n}$${\ textstyle {\ frac {4} {3}} n ^ {3}}$

Invertendo uma matriz

Ao resolver sistemas de equações, b é geralmente tratado como um vector com um comprimento igual à altura da matriz Uma . Na inversão da matriz, no entanto, em vez do vetor b , temos a matriz B , onde B é uma matriz n -by- p , de modo que estamos tentando encontrar uma matriz X (também uma matriz n -by- p ):

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

Podemos usar o mesmo algoritmo apresentado anteriormente para resolver para cada coluna de matriz X . Agora suponha que B seja a matriz de identidade de tamanho n . Segue-se que o resultado X deve ser o inverso de A .

Calculando o determinante

Dada a decomposição LUP de uma matriz quadrada A , o determinante de A pode ser calculado diretamente como ${\ 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} \ direita) \ esquerda (\ prod _ {i = 1} ^ {n} u_ {ii} \ direita).}$

A segunda equação segue do fato de que o determinante de uma matriz triangular é simplesmente o produto de suas entradas diagonais, e que o determinante de uma matriz de permutação é igual a (−1) ^S onde S é o número de trocas de linha na decomposição .

No caso da decomposição LU com pivotamento completo, também é igual ao lado direito da equação acima, se deixarmos S ser o número total de trocas de linha e coluna. ${\ textstyle \ det (A)}$

O mesmo método se aplica prontamente à decomposição LU, definindo P igual à matriz de identidade.

Exemplos de código

Exemplo de código 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;
}

Exemplo de código 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;
    }
}

Exemplo de código 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

Veja também

• Bloco de decomposição LU
• Decomposição de Bruhat
• Decomposição de Cholesky
• Fatoração LU incompleta
• Redução LU
• Decomposição de matriz
• Decomposição QR

Notas

Referências

• Bunch, James R .; Hopcroft, John (1974), "Triangular factorization and inversion by fast matrix multiplation", 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), Introdução a Algoritmos , MIT Press e McGraw-Hill, ISBN 978-0-262-03293-3.
• Golub, Gene H .; Van Loan, Charles F. (1996), Matrix Computations (3ª ed.), 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. Consulte a Seção 3.5. N  - 1
• Householder, Alston S. (1975), The Theory of Matrices in Numerical Analysis , Nova York: Dover Publications , MR  0378371.
• Okunev, Pavel; Johnson, Charles R. (1997), Necessary And Sufficient Conditions For Existence of the LU Factorization of an Arbitrary Matrix , arXiv : math.NA/0506382.
• Poole, David (2006), Linear Algebra: A Modern Introduction (2ª ed.), Canadá: Thomson Brooks / Cole, ISBN 978-0-534-99845-5.
• Pressione, WH; Teukolsky, SA; Vetterling, WT; Flannery, BP (2007), "Section 2.3" , Numerical Recipes: The Art of Scientific Computing (3rd ed.), New York: Cambridge University Press, ISBN 978-0-521-88068-8.
• Trefethen, Lloyd N .; Bau, David (1997), Numerical linear algebra , Philadelphia: Society for Industrial and Applied Mathematics, ISBN 978-0-89871-361-9.

links externos

Referências

• Decomposição LU em MathWorld .
• Decomposição LU no Math-Linux .
• Decomposição LU no Holistic Numerical Methods Institute
• Fatoração da matriz LU . Referência do MATLAB.

Código de computador

• LAPACK é uma coleção de sub-rotinas FORTRAN para resolver problemas densos de álgebra linear
• ALGLIB inclui uma porta parcial do LAPACK para C ++, C #, Delphi, etc.
• Código C ++ , Prof. J. Loomis, University of Dayton
• C code , Mathematics Source Library
• LU em X10

Recursos online

• WebApp resolvendo descritivamente sistemas de equações lineares com Decomposição LU
• Calculadora de matriz , bluebit.gr
• Ferramenta de decomposição LU , uni-bonn.de
• LU Decomposition por Ed Pegg, Jr. , The Wolfram Demonstrations Project , 2007.