Decomposizione LU

In algebra lineare una decomposizione LU, o decomposizione LUP o decomposizione di Doolittle è una fattorizzazione di una matrice in una matrice triangolare inferiore ${\displaystyle L}$, una matrice triangolare superiore ${\displaystyle U}$ e una matrice di permutazione ${\displaystyle P}$. Questa decomposizione è usata in analisi numerica per risolvere un sistema di equazioni lineari, per calcolare l'inversa di una matrice o per calcolare il determinante di una matrice.

Definizione

Sia ${\displaystyle A}$ una matrice invertibile. Allora ${\displaystyle A}$ può essere decomposta come:

${\displaystyle PA=LU}$

dove ${\displaystyle P}$ è una matrice di permutazione, ${\displaystyle L}$ è una matrice triangolare inferiore a diagonale unitaria (${\displaystyle l_{ii}=1,\forall i}$) e ${\displaystyle U}$ è una matrice triangolare superiore.

Idea principale

La decomposizione ${\displaystyle LU}$ è simile all'algoritmo di Gauss. Nell'eliminazione gaussiana si prova a risolvere l'equazione matriciale:

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

Il processo di eliminazione produce una matrice triangolare superiore ${\displaystyle U}$ e trasforma il vettore ${\displaystyle b}$ in ${\displaystyle b'}$

${\displaystyle Ux=b'}$

Poiché ${\displaystyle U}$ è una matrice triangolare superiore, questo sistema di equazioni si può risolvere facilmente tramite sostituzione all'indietro.

Durante la decomposizione ${\displaystyle LU}$, però, ${\displaystyle b}$ non è trasformato e l'equazione può essere scritta come:

${\displaystyle Ax=LUx=b}$

così si può riutilizzare la decomposizione se si vuole risolvere lo stesso sistema per un differente ${\displaystyle b}$.

Nel caso più generale, nel quale la fattorizzazione della matrice comprende anche l'utilizzo di scambi di riga nella matrice, viene introdotta anche una matrice di permutazione ${\displaystyle P}$, ed il sistema diventa:

${\displaystyle {\begin{cases}Ly=Pb\\Ux=y\end{cases}}}$

La risoluzione di questo sistema permette la determinazione del vettore ${\displaystyle x}$ cercato.

Algoritmo

Applicando delle serie di trasformazioni elementari di matrice (cioè moltiplicazioni di ${\displaystyle A}$ a sinistra) si costruisce una matrice triangolare superiore ${\displaystyle U}$ che parte da ${\displaystyle A}$. Questo metodo è chiamato metodo di Gauss. Queste trasformazioni elementari di matrice sono tutte delle trasformazioni lineari di tipo combinatorio (il terzo tipo elencato nella voce "Matrice elementare"). Si supponga che ${\displaystyle T}$ sia il prodotto di N trasformazioni ${\displaystyle T_{N}\dots T_{2}T_{1}=T}$, allora la matrice triangolare superiore è:

${\displaystyle TA=T_{N}\dots T_{2}T_{1}A=\colon U}$

L'inversa della matrice ${\displaystyle T}$ è:

${\displaystyle T^{-1}=T_{1}^{-1}T_{2}^{-1}\dots T_{N}^{-1}}$

Come l'algoritmo di Gauss usa solo la terza forma dei tre tipi di trasformazioni elementari di matrice rendendo ${\displaystyle A}$ triangolare superiore, si può dedurre che tutte le ${\displaystyle T_{i}^{-1}}$ sono triangolari inferiori (vedi trasformazioni elementari di matrice). Essendo un prodotto di ${\displaystyle T_{i}^{-1}}$ anche:

${\displaystyle T^{-1}=T_{1}^{-1}T_{2}^{-1}\dots T_{N}^{-1}=\colon L}$

è triangolare inferiore.

Si ha quindi la decomposizione della matrice ${\displaystyle A}$ nel prodotto di ${\displaystyle L}$ e ${\displaystyle U}$:

${\displaystyle LU=T^{-1}TA=A}$

Applicazioni

Matrici inverse

La fattorizzazione ${\displaystyle PA=LU}$ viene anche usata per calcolare la matrice inversa di ${\displaystyle A}$. Infatti:

${\displaystyle PA=LU\Rightarrow A=P^{-1}LU}$

da cui:

${\displaystyle A^{-1}=(P^{-1}LU)^{-1}=U^{-1}L^{-1}P}$

Determinante

Un altro utilizzo di questa decomposizione è per il calcolo del determinante della matrice ${\displaystyle A}$. Infatti:

${\displaystyle A=P^{-1}LU\Rightarrow \det(A)=\det(P^{-1}LU)}$

quindi per il teorema di Binet:

${\displaystyle \det(A)=\det(P^{-1})\det(L)\det(U)}$

Sapendo che il determinante di una matrice di permutazione vale ${\displaystyle 1}$ se questa corrisponde ad un numero pari di permutazioni, ${\displaystyle -1}$ altrimenti, e che ${\displaystyle \det A^{-1}={\frac {1}{\det A}}}$, si ottiene che:

${\displaystyle \det(A)=(-1)^{S}\det(L)\det(U)}$

Sapendo poi che il determinante di una matrice triangolare è dato dal prodotto dei termini sulla sua diagonale principale, si ha che:

${\displaystyle \det(A)=\left(-1\right)^{S}\prod _{i=1}^{n}{u_{ii}l_{ii}}}$

ma sapendo anche che i termini sulla diagonale principale di ${\displaystyle L}$ sono tutti ${\displaystyle 1}$, quindi si può infine concludere:

${\displaystyle \det(A)=\left(-1\right)^{S}\prod _{i=1}^{n}u_{ii}}$

dove ${\displaystyle S}$ indica il numero di scambi di riga effettuati nel processo (implicitamente indicati nella matrice ${\displaystyle P}$) ed i termini ${\displaystyle u_{ij}}$ e ${\displaystyle l_{ij}}$ indicano il termine in riga ${\displaystyle i}$ e colonna ${\displaystyle j}$ rispettivamente delle matrici ${\displaystyle U}$ e ${\displaystyle L}$.

Codice in C

/* INPUT: A - vettore di puntatori alle righe della matrice quadrata di dimensione N
 *        Tol - Callore della tolleranza minima per identificare quando la matrice è prossima a degenerare
 * OUTPUT: La matrice A è cambiata, contiene sia le matrici L-E e U tali che A = (L-E)+U e P*A = L*U
 *         La matrice di permutazione non è salvata in memoria come una matrice, ma in un vettore
           di interi di dimensione N+1   
 *         contenente gli indici delle colonne in cui la matrice ha come elementi "1". L'ultimo elemento P[N]=S+N, 
 *         dove S è il numero delle righe scambiate necessarie per il calcolo del determinante, 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; //Matrice di permutazione unaria, P[N] inizializzato con 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; //La matrice è degenerata

        if (imax != i) {
            //pivoting P
            j = P[i];
            P[i] = P[imax];
            P[imax] = j;

            //pivoting delle righe di A
            ptr = A[i];
            A[i] = A[imax];
            A[imax] = ptr;

            //Conteggio dei pivot partendo da N per il calcolo del determinante
            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;  //Decomposizione conclusa
}

/* INPUT: A,P matrici riempite nella funzione LUPDecompose; b - vettore; N - dimensione
 * OUTPUT: x - vettore soluzione di 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] = x[i] / A[i][i];
    }
}

/* INPUT: A,P matrici riempite nella funzione LUPDecompose; N - dimensione
 * OUTPUT: IA è l'inversa della matrice iniziale
 */
void LUPInvert(double **A, int *P, int N, double **IA) {
  
    for (int j = 0; j < N; j++) {
        for (int i = 0; i < N; i++) {
            if (P[i] == j) 
                IA[i][j] = 1.0;
            else
                IA[i][j] = 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] = IA[i][j] / A[i][i];
        }
    }
}

/* INPUT: A,P matrici riempite nella funzione LUPDecompose; N - dimensione. 
 * OUTPUT: La funzione restituisce il valore del determinante della matrice iniziale
 */
double LUPDeterminant(double **A, int *P, int N) {

    double det = A[0][0];

    for (int i = 1; i < N; i++)
        det *= A[i][i];

    if ((P[N] - N) % 2 == 0)
        return det; 
    else
        return -det;
}

Codice in 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
            // Calcolo delle soluzioni di 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;
            }
            // Calcolo delle soluzioni di 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;
        }
}

Codice in Matlab

function x = SolveLinearSystem(A, b)
    n = length(b);
    x = zeros(n, 1);
    y = zeros(n, 1);
    % Decomposizione della matrice per mezzo del metodo Doolittle
    for i = 1:1:n
        for j = 1:1:(i - 1)
            alpha = A(i,j);
            for k = 1:1:(j - 1)
                alpha = alpha - A(i,k)*A(k,j);
            end
            A(i,j) = alpha/A(j,j);
        end
        for j = i:1:n
            alpha = A(i,j);
            for k = 1:1:(i - 1)
                alpha = alpha - A(i,k)*A(k,j);
            end
            A(i,j) = alpha;
        end
    end
    % A = L+U-I
    % calcolo delle soluzioni di Ly = b
    for i = 1:1:n
        alpha = 0;
        for k = 1:1:i
            alpha = alpha + A(i,k)*y(k);
        end
        y(i) = b(i) - alpha;
    end
    % calcolo delle soluzioni di Ux = y
    for i = n:(-1):1
        alpha = 0;
        for k = (i + 1):1:n
            alpha = alpha + A(i,k)*x(k);
        end
        x(i) = (y(i) - alpha)/A(i, i);
    end    
end

Bibliografia

• (EN) Bunch, James R.; Hopcroft, John (1974), "Triangular factorization and inversion by fast matrix multiplication", Mathematics of Computation 28: 231–236, ISSN 0025-5718.
• (EN) Cormen, Thomas H.; Leiserson, Charles E.; Rivest, Ronald L.; Stein, Clifford (2001), Introduction to Algorithms, MIT Press and McGraw-Hill, ISBN 978-0-262-03293-3.
• (EN) Golub, Gene H.; Van Loan, Charles F. (1996), Matrix Computations (3rd ed.), Baltimore: Johns Hopkins, ISBN 978-0-8018-5414-9.
• (EN) Horn, Roger A.; Johnson, Charles R. (1985), Matrix Analysis, Cambridge University Press, ISBN 0-521-38632-2. See Section 3.5.
• (EN) Householder, Alston S. (1975), The Theory of Matrices in Numerical Analysis, New York: Dover Publications, MR 0378371.
• (EN) Okunev, Pavel; Johnson, Charles R. (1997), Necessary And Sufficient Conditions For Existence of the LU Factorization of an Arbitrary Matrix, arXiv:math.NA/0506382.
• (EN) Poole, David (2006), Linear Algebra: A Modern Introduction (2nd ed.), Canada: Thomson Brooks/Cole, ISBN 0-534-99845-3.
• (EN) Press, 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.
• (EN) Trefethen, Lloyd N.; Bau, David (1997), Numerical linear algebra, Philadelphia: Society for Industrial and Applied Mathematics, ISBN 978-0-89871-361-9.

Voci correlate

• Analisi numerica
• Decomposizione di Cholesky
• Decomposizione di Crout
• Decomposizione di una matrice
• Decomposizione QR
• Matrice triangolare

Altri progetti

Altri progetti

• Wikimedia Commons

• Wikimedia Commons contiene immagini o altri file su Decomposizione LU

Collegamenti esterni

• (EN) Eric W. Weisstein, Decomposizione LU, su MathWorld, Wolfram Research.
• (EN) LU decomposition on Math-Linux.
• (EN) Module for LU Factorization with Pivoting, Prof. J. H. Mathews, California State University, Fullerton
• (EN) LU decomposition at Holistic Numerical Methods Institute
• (EN) LU matrix factorization. MATLAB reference.
• (EN) WebApp descriptively solving systems of linear equations with LU Decomposition, su sole.ooz.ie. URL consultato il 10 aprile 2014 (archiviato dall'url originale il 25 aprile 2011).
• (EN) Matrix Calculator, bluebit.gr
• (EN) LU Decomposition Tool, uni-bonn.de
• (EN) LU Decomposition by Ed Pegg, Jr., The Wolfram Demonstrations Project, 2007.