Fattorizzazione lU (eliminazione gaussiana)

Eliminazione gaussiana

Considero la matrice $A$$A$ e faccio operazioni sulle sue righe.

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

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

Risolvo direttamente l'eliminazione di Gauss sulla matrice che ha come colonne le colonne di $A$$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$$A=LU$, devo risolvere i sistemi $L{\boldsymbol {y}}={\boldsymbol {b}}$$L{\boldsymbol {y}}={\boldsymbol {b}}$ e poi $U{\boldsymbol {x}}={\boldsymbol {y}}$$U{\boldsymbol {x}}={\boldsymbol {y}}$, spezzo 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}}$$A \boldsymbol{x} = \boldsymbol{b}$ in un sistema della forma $U{\boldsymbol {x}}={\boldsymbol {c}}$$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}$$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$$k+1$-esima all'$n$$n$-esima.

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

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

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

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

e sottraggo ad $A_{i}$$A_i$ la riga $M_{ik}*A_{k}$$M_{ik}*A_{k}$, in modo che

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

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

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

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

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

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

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

Matrice elementare di Gauss

Definizione 2.11

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

1. Nel blocco $k-1\times k-1$$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}}$$

   $${\begin{array}{c}1\\-M_{k+1,k}\\-M_{k+2,k}\\\dots \\-M_{n,k}\end{array}}$$
   dove gli $M_{ik}$$M_{ik}$ sono i coefficienti definiti prima, mentre la prima riga di questo blocco è il vettore ${\boldsymbol {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}}$$

$${\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}$$M_{k}=I-{\boldsymbol {w}}_{k}*{\boldsymbol {e}}_{k}^{t}$, dove ${\boldsymbol {w}}_{k}$${\boldsymbol {w}}_{k}$ è un vettore colonna con le prime k componenti uguali a 0, e le restanti uguali a $M_{ik}$$M_{ik}$ e ${\boldsymbol {e}}_{k}$${\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^{(1)}=A$, $A^{(k+1)}=M_{k}A^{(k)}$$A^{(k+1)}=M_{k}A^{(k)}$ per $k=1,\dots ,n-1$$k=1,\dots ,n-1$, dove $A^{(k)}$$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$$U$ la matrice triangolare superiore ottenuta all'ultimo passo, si ha:

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

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

Per poter scrivere $A=LU$$A=LU$, ci si chiede se il prodotto delle $M_{i}$$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$$

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

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

Osservazione 2.3

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

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

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

Allora

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

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

$$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}$${\boldsymbol {e}}_{k}^{t}{\boldsymbol {w}}_{k}$ è il vettore nullo, gli altri due termini sono opposti e rimane $M_{k}^{-1}M_{k}=I$$M_{k}^{-1}M_{k}=I$.

Considero allora $L=M_{1}^{-1}M_{2}^{-1}\dots M_{n-1}^{-1}$$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+{\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}$$

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

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

e ${\boldsymbol {e}}_{1}^{t}{\boldsymbol {w}}_{2}=0$${\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$$

$${\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}$$

$$(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}$${\boldsymbol {w}}_{k}{\boldsymbol {e}}_{k}^{t}$ si sommano una accanto all'altra, e si ottiene la matrice triangolare inferiore $L$$L$ cercata, con entrate uguali a 1 sulla diagonale, zeri nella parte superiore, e i moltiplicatori nella parte inferiore. (i moltiplicatori in $L$$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}}$$L{\boldsymbol {y}}={\boldsymbol {b}}$. Eseguo, per $k=1,\dots ,n$$k=1,\dots ,n$ e per $i=k+1,\dots ,n$$i=k+1,\dots ,n$:

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

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

$$y_{1}=b_{1}$$

$$y_{1}=b_{1}$$

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

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

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

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

$$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}$$n^{3}$ operazioni, mentre la risoluzione vera e propria ne richiede $n^{2}$$n^2$.

Riepilogo

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

1. definiamo i moltiplicatori
   $$M_{ik}={\frac {|A_{ik}|}{|A_{kk}|}}$$

   $$M_{ik}={\frac {|A_{ik}|}{|A_{kk}|}}$$
2. per $j=k+1,\dots ,n$$j=k+1,\dots ,n$
   $$a_{ij}^{(k+1)}=a_{ij}^{(k)}-M_{ik}a_{kj}^{(k)}$$

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

   $$b_{i}^{(k+1)}=b_{i}^{(k)}-M_{ik}b_{k}^{(k)}$$
4. passo quindi dal sistema $A{\boldsymbol {x}}={\boldsymbol {b}}$$A \boldsymbol{x} = \boldsymbol{b}$ al sistema $U{\boldsymbol {x}}={\boldsymbol {c}}$$U{\boldsymbol {x}}={\boldsymbol {c}}$ con ${\boldsymbol {c}}$${\boldsymbol {c}}$ soluzione di $L{\boldsymbol {y}}={\boldsymbol {b}}$$L{\boldsymbol {y}}={\boldsymbol {b}}$.
5. $U$$U$ è la matrice ottenuta dopo l'eliminazione gaussiana, mentre $L$$L$ ha entrate uguali a 1 sulla diagonale principale e ha i coefficienti $M_{ik}$$M_{ik}$ nella parte inferiore.

Problema dell'esistenza della fattorizzazione

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

Teorema 2.5

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

 

Nel momento del calcolo del moltiplicatore $M_{ik}$$M_{ik}$, in aritmetica esatta si hanno problemi se $A_{kk}=0$$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}$$A \in \mathbb C_{n,n}$, allora esiste $P$$P$ matrice di permutazione tale che $PA$$PA$ è fattorizzabile nel prodotto $LU$$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$$A_{ii}=A_{kk}=1$ e $A_{ik}=A_{ki}=0$$A_{ik}=A_{ki}=0$, si inverte il loro valore, ponendo $A_{ii}=A_{kk}=0$$A_{ii}=A_{kk}=0$ e $A_{ik}=A_{ki}=1$$A_{ik}=A_{ki}=1$.