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 , una matrice triangolare superiore e una matrice di permutazione . 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.

DefinizioneModifica

Sia   una matrice invertibile. Allora   può essere decomposta come:

 

dove   è una matrice di permutazione,   è una matrice triangolare inferiore a diagonale unitaria ( ) e   è una matrice triangolare superiore.

Idea principaleModifica

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

 

Il processo di eliminazione produce una matrice triangolare superiore   e trasforma il vettore   in  

 

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

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

 

così si può riutilizzare la decomposizione se si vuole risolvere lo stesso sistema per un differente  .

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  , ed il sistema diventa:

 

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

AlgoritmoModifica

Applicando delle serie di trasformazioni elementari di matrice (cioè moltiplicazioni di   a sinistra) si costruisce una matrice triangolare superiore   che parte da  . 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   sia il prodotto di N trasformazioni  , allora la matrice triangolare superiore è:

 

L'inversa della matrice   è:

 

Come l'algoritmo di Gauss usa solo la terza forma dei tre tipi di trasformazioni elementari di matrice rendendo   triangolare superiore, si può dedurre che tutte le   sono triangolari inferiori (vedi trasformazioni elementari di matrice). Essendo un prodotto di   anche:

 

è triangolare inferiore.

Si ha quindi la decomposizione della matrice   nel prodotto di   e  :

 

ApplicazioniModifica

Matrici inverseModifica

La fattorizzazione   viene anche usata per calcolare la matrice inversa di  . Infatti:

 

da cui:

 

DeterminanteModifica

Un altro utilizzo di questa decomposizione è per il calcolo del determinante della matrice  . Infatti:

 

quindi per il teorema di Binet:

 

Sapendo che il determinante di una matrice di permutazione vale   se questa corrisponde ad un numero pari di permutazioni,   altrimenti, e che  , si ottiene che:

 

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

 

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

 

dove   indica il numero di scambi di riga effettuati nel processo (implicitamente indicati nella matrice  ) ed i termini   e   indicano il termine in riga   e colonna   rispettivamente delle matrici   e  .

Codice in CModifica

/* 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#Modifica

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 MatlabModifica

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

BibliografiaModifica

  • (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 correlateModifica

Collegamenti esterniModifica

  Portale Matematica: accedi alle voci di Wikipedia che trattano di matematica