Calcolo numerico/Fattorizzazione LU (eliminazione gaussiana)

Eliminazione gaussiana

Consideriamo la matrice $A$ e facciamo operazioni sulle sue righe.

$A{\boldsymbol {x}}=[{\boldsymbol {b}}]^{i},\quad i=1,\dots ,n$

Risolviamo direttamente l'eliminazione di Gauss sulla matrice che ha come colonne le colonne di $A$ e poi i termini noti. Se i termini non sono noti a priori possono esserci dei problemi, e si può applicare il metodo della potenza inversa, che serve per calcolare l'autovalore di modulo minimo.

Dopo aver fattorizzato $A=LU$, bisogna risolvere i sistemi $L{\boldsymbol {y}}={\boldsymbol {b}}$ e poi $U{\boldsymbol {x}}={\boldsymbol {y}}$. Spezziamo quindi il problema in due parti, separando la fattorizzazione della matrice dalla trattazione vera e propria della soluzione.

Descrizione dell'algoritmo

L'eliminazione gaussiana trasforma un sistema della forma $A{\boldsymbol {x}}={\boldsymbol {b}}$ in un sistema della forma $U{\boldsymbol {x}}={\boldsymbol {c}}$.

Supponiamo di essere al $k$-esimo passo di fattorizzazione, cioè supponiamo che la parte della matrice superiore alla riga $A_{k}$ sia già stata ridotta a una triangolare superiore. A ogni passo, le operazioni sulle righe della matrice modificano anche il blocco quadrato in basso a destra, con le righe dalla $k+1$-esima all'$n$-esima.

Sottraiamo alla $i$-esima equazione la $k$-esima moltiplicata per un opportuno coefficiente.

Descrivendo l'algoritmo, per $k=1,\dots ,n-1$ (passi di fattorizzazione), per $i=k+1,\dots ,n$ (righe successive alla $k$-esima), definisco l'elemento

$M_{ik}={\frac {A_{ik}^{(k)}}{A_{kk}^{(k)}}}$

e sottraggo a $M_{ik}*A_{k}$, in modo che

$A_{ik}^{(k+1)}=0$

Ad esempio, nel passo $k=1$, prendiamo la prima riga e la sottraiamo a tutte le altre moltiplicata per il rispettivo coefficiente $M_{ik}$. Nel passo successivo lasciamo invariate le righe precedenti alla $k$-esima, e così via.

I coefficienti nel blocco quadrato degli $A_{ij}$ con $k+1<j<n,k+1<i<n$ in basso a destra, che vengono modificati a ogni passo $k$, diventano:

$A_{ij}^{(k+1)}=A_{ij}^{(k)}-M_{ik}*A_{kj}^{(k)}$

Analogamente per il termine noto:

$b_{i}^{(k+1)}=b_{i}^{(k)}-M_{ik}b_{k}^{(k)}$

Matrice elementare di Gauss

Definizione

Sia $M_{k}$ una matrice triangolare inferiore, con $\det M_{k}=1$, divisa in quattro blocchi e tale che:

1. nel blocco $k-1\times k-1$ in alto a sinistra c'è l'identità,
2. i blocchi di nord-est e di sud-ovest sono nulli,
3. la matrice in basso a destra ha la prima sottocolonna della forma

${\begin{array}{c}1\\-M_{k+1,k}\\-M_{k+2,k}\\\dots \\-M_{n,k}\end{array}}$

dove gli $M_{ik}$ sono i coefficienti definiti prima, mentre la prima riga di questo blocco è il vettore e 1 ${\boldsymbol {e}}_{1}$, e nel quadrato in basso c'è l'identità.

In forma sintetica:

${\begin{array}{|ccc|c|cc|}\hline \\&&&&&\\&I_{k-1}&&&0&\\&&&&&\\\hline \\&&&1&0&0\\\hline {5:6}\\&0&&-M_{k+1,k}&I_{k}&\\&&&-M_{n,k}&&\\\hline \end{array}}$

In maniera compatta: $M_{k}=I-{\boldsymbol {w}}_{k}*{\boldsymbol {e}}_{k}^{t}$, dove ${\boldsymbol {w}}_{k}$ è un vettore colonna con le prime $k$ componenti uguali a 0, e le restanti uguali a $M_{ik}$ e ${\boldsymbol {e}}_{k}$ è il $k$-esimo vettore della base canonica. Il prodotto tra il vettore colonna e il vettore riga è una matrice.

Rappresentazione compatta dell'eliminazione di Gauss

Poniamo $A^{(1)}=A$, $A^{(k+1)}=M_{k}A^{(k)}$ per $k=1,\dots ,n-1$, dove $A^{(k)}$ indica la matrice ottenuta al passo $k$-esimo. Ogni passo dell'eliminazione di Gauss può essere quindi rappresentato in modo compatto come un prodotto tra due matrici. Allora, chiamando $U$ la matrice triangolare superiore ottenuta all'ultimo passo, si ha:

$U=M_{n-1}\dots M_{2}M_{1}A$

Per poter scrivere $A=LU$, ci si chiede se il prodotto delle $M_{i}$ è invertibile, e questo è vero perché è prodotto di matrici invertibili. Allora si può scrivere:

$A=(M_{n-1}M_{n-2}\dots M_{2}M_{1})^{-1}U$

ed esplicitando la matrice inversa del prodotto:

$A=M_{1}^{-1}M_{2}^{-1}\dots M_{n-1}^{-1}U$

Osservazione 2.3

Siccome $L$ è la matrice triangolare inferiore con entrate uguali a 1 sulla diagonale, il suo determinante, prodotto delle entrate diagonali, è 1, e quindi

$\det A=\det(LU)=\det L*\det U=\det U$

allora la fattorizzazione della matrice può essere considerata anche come un modo per calcolare il determinante.

Teniamo conto del fatto che:

$M_{k}=I-{\boldsymbol {w}}_{k}{\boldsymbol {e}}_{k}^{t}$

Allora

$M_{k}^{-1}=I+{\boldsymbol {w}}_{k}{\boldsymbol {e}}_{k}^{t}$

infatti si verifica che

$M_{k}^{-1}M_{k}=(I+{\boldsymbol {w}}_{k}{\boldsymbol {e}}_{k}^{t})(I-{\boldsymbol {w}}_{k}{\boldsymbol {e}}_{k}^{t})=I+{\boldsymbol {w}}_{k}{\boldsymbol {e}}_{k}^{t}-{\boldsymbol {w}}_{k}{\boldsymbol {e}}_{k}^{t}-{\boldsymbol {w}}_{k}{\boldsymbol {e}}_{k}^{t}{\boldsymbol {w}}_{k}{\boldsymbol {e}}_{k}^{t}$

ma ${\boldsymbol {e}}_{k}^{t}{\boldsymbol {w}}_{k}$ è il vettore nullo, gli altri due termini sono opposti e rimane $M_{k}^{-1}M_{k}=I$.

Consideriamo allora $L=M_{1}^{-1}M_{2}^{-1}\dots M_{n-1}^{-1}$, che si riscrive come

$=(I+{\boldsymbol {w}}_{1}{\boldsymbol {e}}_{1}^{t})*(I+{\boldsymbol {w}}_{2}{\boldsymbol {e}}_{2}^{t})\dots (I+{\boldsymbol {w}}_{k-1}{\boldsymbol {e}}_{k-1}^{t})$

$=I+\sum _{i=1}^{n-1}{\boldsymbol {w}}_{i}{\boldsymbol {e}}_{i}^{t}$

L'ultimo passaggio vale perché, ad esempio

$(I+{\boldsymbol {w}}_{1}{\boldsymbol {e}}_{1}^{t})(I+{\boldsymbol {w}}_{2}{\boldsymbol {e}}_{2}^{t})=I+{\boldsymbol {w}}_{1}{\boldsymbol {e}}_{1}^{t}+{\boldsymbol {w}}_{2}{\boldsymbol {e}}_{2}^{t}+{\boldsymbol {w}}_{1}{\boldsymbol {e}}_{1}^{t}{\boldsymbol {w}}_{2}{\boldsymbol {e}}_{2}^{t}$

${\boldsymbol {e}}_{1}^{t}{\boldsymbol {w}}_{2}=0$ come prima, perché hanno zeri in posizioni complementari. In generale

${\boldsymbol {e}}_{i}^{t}{\boldsymbol {w}}_{k}=0\forall i\leq k$

i termini misti si annullano, e quindi rimane

$(I+{\boldsymbol {w}}_{1}{\boldsymbol {e}}_{1}^{t})(I+{\boldsymbol {w}}_{2}{\boldsymbol {e}}_{2}^{t})=I+{\boldsymbol {w}}_{1}{\boldsymbol {e}}_{1}^{t}+{\boldsymbol {w}}_{2}{\boldsymbol {e}}_{2}^{t}$

e, procedendo così fino al passo $n$-esimo, le colonne ${\boldsymbol {w}}_{k}{\boldsymbol {e}}_{k}^{t}$ si sommano una accanto all'altra, e si ottiene la matrice triangolare inferiore $L$ cercata, con entrate uguali a 1 sulla diagonale, zeri nella parte superiore, e i moltiplicatori nella parte inferiore (i moltiplicatori in $L$ sono positivi).

Operazioni sul termine noto

Le operazioni che l'eliminazione gaussiana fa sul termine noto equivalgono alla risoluzione del sistema $L{\boldsymbol {y}}={\boldsymbol {b}}$. Eseguo, per $k=1,\dots ,n$ e per $i=k+1,\dots ,n$:

$y_{i}=b_{i}-L_{ik}y_{k}=b_{i}-M_{ik}y_{k}.$

Esempio 2.6

Ad esempio, nel caso di una matrice $4\times 4$,

$y_{1}=b_{1}$

$y_{2}=b_{2}-L_{21}y_{1}$

$y_{3}=b_{3}-L_{31}y_{1}-L_{32}y_{2}$

$y_{4}=b_{4}-L_{41}y_{1}-L_{42}y_{2}-L_{43}y_{3}$

Al primo passo di eliminazione gaussiana, il primo coefficiente è già definitivo mentre tutti gli altri vengono modificati. Al secondo passo anche il secondo coefficiente risulta definitivo, e così via. Al terzo passo, si lavora solo sull'ultima componente.

Si ha un aggiornamento progressivo per colonne.

La fattorizzazione ha un costo computazionale di $n^{3}$ operazioni, mentre la risoluzione vera e propria ne richiede $n^{2}$.

Riepilogo

Vogliamo risolvere il sistema $A{\boldsymbol {x}}={\boldsymbol {b}}$, e nell'algoritmo dell'eliminazione gaussiana, per $k=1,\dots ,n-1$, per $i=k+1,\dots ,n$, compiamo le seguenti operazioni:

1. definiamo i moltiplicatori

$M_{ik}={\frac {|A_{ik}|}{|A_{kk}|}}$

2. per $j=k+1,\dots ,n$

$a_{ij}^{(k+1)}=a_{ij}^{(k)}-M_{ik}a_{kj}^{(k)}$

3. sul termine noto

$b_{i}^{(k+1)}=b_{i}^{(k)}-M_{ik}b_{k}^{(k)}$

4. passo quindi dal sistema $A{\boldsymbol {x}}={\boldsymbol {b}}$ al sistema $U{\boldsymbol {x}}={\boldsymbol {c}}$ con ${\boldsymbol {c}}$ soluzione di $L{\boldsymbol {y}}={\boldsymbol {b}}$.

5. $U$ è la matrice ottenuta dopo l'eliminazione gaussiana, mentre $L$ ha entrate uguali a 1 sulla diagonale principale e ha i coefficienti $M_{ik}$ nella parte inferiore.

Problema dell'esistenza della fattorizzazione

Data una matrice quadrata $A$, chiamiamo sottomatrici principali di $A$ le matrici $A_{k}$, ottenute privando $A$ della $k$-esima riga e della $k$-esima colonna. Si dicono matrici principali di testa le sottomatrici che contengono il coefficiente $A_{11}$.

Teorema 2.5

Data una matrice quadrata $A$ in $\mathbb {C} _{n,n}$, se le matrici $A_{k}$ sono non singolari per $k=1\dots ,n$, allora esiste unica la fattorizzazione $A=LU$.

Nel momento del calcolo del moltiplicatore $M_{ik}$, in aritmetica esatta si hanno problemi se $A_{kk}=0$. Siccome la matrice è non singolare, esisterà una riga in cui il termine in posizione $k$ è diverso da 0, allora, per evitare problemi, si scambia questa riga con la $k$-esima.

Teorema 2.6

Data $A\in \mathbb {C} _{n,n}$, allora esiste $P$ matrice di permutazione tale che $PA$ è fattorizzabile nel prodotto $LU$.

Una matrice di permutazione che scambia le righe $i$-esima e $k$-esima si ottiene scambiando tali righe nella matrice identica, in particolare, considerando i quattro coefficienti $A_{ii}=A_{kk}=1$ e $A_{ik}=A_{ki}=0$, si inverte il loro valore, ponendo $A_{ii}=A_{kk}=0$ e $A_{ik}=A_{ki}=1$.