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 $L$ , una matrice triangolare superiore $U$ e una matrice di permutazione $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 $A$ una matrice invertibile. Allora $A$ può essere decomposta come:

PA=LU

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

Idea principale

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

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

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

Ux=b'

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

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

Ax=LUx=b

così si può riutilizzare la decomposizione se si vuole risolvere lo stesso sistema per un differente $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 $P$ , ed il sistema diventa:

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

La risoluzione di questo sistema permette la determinazione del vettore $x$ cercato.

Algoritmo

Applicando delle serie di trasformazioni elementari di matrice (cioè moltiplicazioni di $A$ a sinistra) si costruisce una matrice triangolare superiore $U$ che parte da $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 $T$ sia il prodotto di N trasformazioni $T_{N}\dots T_{2}T_{1}=T$ , allora la matrice triangolare superiore è:

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

L'inversa della matrice $T$ è:

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 $A$ triangolare superiore, si può dedurre che tutte le $T_{i}^{-1}$ sono triangolari inferiori (vedi trasformazioni elementari di matrice). Essendo un prodotto di $T_{i}^{-1}$ anche:

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

è triangolare inferiore.

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

LU=T^{-1}TA=A

Applicazioni

Matrici inverse

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

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

da cui:

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 $A$ . Infatti:

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

quindi per il teorema di Binet:

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

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

$\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:

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

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

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

dove $S$ indica il numero di scambi di riga effettuati nel processo (implicitamente indicati nella matrice $P$ ) ed i termini $u_{ij}$ e $l_{ij}$ indicano il termine in riga $i$ e colonna $j$ rispettivamente delle matrici $U$ e $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

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.