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.

Idea intuitiva

Si immagini di dover risolvere un sistema di equazioni lineari come un puzzle complesso. La decomposizione LU è come scomporre questo puzzle in due puzzle più semplici da risolvere in sequenza.

Una matrice $A$ può essere "spezzata" in due parti più managevoli:

• una matrice triangolare inferiore $L$ (Lower) - ha tutti zeri sopra la diagonale;
• una matrice triangolare superiore $U$ (Upper) - ha tutti zeri sotto la diagonale.

È come dire: invece di risolvere direttamente $Ax=b$, risolviamo prima $Ly=b$ (facile, perché $L$ è triangolare), poi $Ux=y$ (anche questo facile, perché $U$ è triangolare).

Perché è utile?

La decomposizione LU è vantaggiosa per tre motivi principali:

Riutilizzo computazionale: se si deve risolvere lo stesso sistema con diversi termini noti $b$, si fa la decomposizione una volta sola e poi si riutilizza. È come assemblare una volta per tutte gli attrezzi giusti, poi usarli per molti lavori simili.

Efficienza numerica: risolvere due sistemi triangolari è molto più veloce che risolvere direttamente il sistema originale. Le matrici triangolari si risolvono "a cascata" dall'alto verso il basso o viceversa.

Stabilità: con il pivoting (scambi di righe), il metodo diventa numericamente stabile anche per matrici "difficili".

Come funziona l'algoritmo (versione semplificata)

L'algoritmo è sostanzialmente l'eliminazione gaussiana, ma invece di dimenticare i passaggi intermedi, vengono inlcusi come parte dell'output:

Passo 1: Trasformare la matrice $A$ in triangolare superiore $U$ usando l'eliminazione gaussiana.

Passo 2: Includere tutti i "moltiplicatori" usati nell'eliminazione in una matrice triangolare inferiore $L$.

Passo 3: Se servono scambi di righe per la stabilità numerica, tenerne traccia con una matrice di permutazione $P$.

Risultato: $PA=LU$, dove $P$ rappresenta gli scambi di righe.

Esempio

Si consideri il sistema:

${\begin{cases}2x+3y=8\\4x+7y=18\end{cases}}$

In forma matriciale:

${\begin{pmatrix}2&3\\4&7\end{pmatrix}}{\begin{pmatrix}x\\y\end{pmatrix}}={\begin{pmatrix}8\\18\end{pmatrix}}.$

La decomposizione LU produrrebbe:

• $L={\begin{pmatrix}1&0\\2&1\end{pmatrix}}$ (il moltiplicatore 2 viene "ricordato");
• $U={\begin{pmatrix}2&3\\0&1\end{pmatrix}}$ (matrice triangolare superiore).

Applicazioni pratiche

La decomposizione LU non è solo teoria accademica, ma ha applicazioni concrete:

Ingegneria strutturale: analisi di telai e travi, dove lo stesso sistema strutturale deve essere analizzato sotto carichi diversi.

Circuiti elettrici: simulazione di reti elettriche complesse, dove la topologia del circuito resta fissa ma cambiano i valori di corrente e tensione.

Computer grafica: trasformazioni geometriche e rendering 3D, dove le stesse matrici di trasformazione vengono applicate a molti punti.

Analisi finanziaria: modelli di portafoglio e pricing di derivati, dove la struttura di correlazione tra asset rimane stabile.

Simulazioni fisiche: fluidodinamica, meccanica dei solidi, dove le equazioni differenziali vengono discretizzate in grandi sistemi lineari che cambiano nel tempo ma mantengono una struttura simile.

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,$ per ogni $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:

$Ax=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$

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

implica

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

dove $s$ indica il numero di scambi di riga effettuati nel processo (implicitamente indicati nella matrice $P$). 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}},$

dove i termini $u_{ij}$ e $l_{ij}$ indicano l'elemento nella riga $i$ e colonna $j$ rispettivamente delle matrici $U$ e $L.$ 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}.$

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, LU Decomposition, 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.

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, LU Decomposition, 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.