Calcolo numerico/Fattorizzazione LU (eliminazione gaussiana)

Eliminazione gaussiana

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

$A𝒙=[𝒃{]}^i,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𝒚=𝒃$ e poi $U𝒙=𝒚$. 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𝒙=𝒃$ in un sistema della forma $U𝒙=𝒄$.

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}\ast 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}\ast A_{kj}^{(k)}$

Analogamente per il termine noto:

$b_i^{(k+1)}=b_i^{(k)}-M_{ik}{b}_k^{(k)}$

Matrice elementare di Gauss

In forma sintetica:

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

In maniera compatta: $M_k=I-{𝒘}_k\ast {𝒆}_k^t$, dove ${𝒘}_k$ è un vettore colonna con le prime $k$ componenti uguali a 0, e le restanti uguali a $M_{ik}$ 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

$detA=det(LU)=detL\ast detU=detU$

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

Teniamo conto del fatto che:

$M_k=I-{𝒘}_k{𝒆}_k^t$

Allora

$M_k^{-1}=I+{𝒘}_k{𝒆}_k^t$

infatti si verifica che

$M_k^{-1}{M}_k=(I+{𝒘}_k{𝒆}_k^t)(I-{𝒘}_k{𝒆}_k^t)=I+{𝒘}_k{𝒆}_k^t-{𝒘}_k{𝒆}_k^t-{𝒘}_k{𝒆}_k^t{𝒘}_k{𝒆}_k^t$

ma ${𝒆}_k^t{𝒘}_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+{𝒘}_1{𝒆}_1^t)\ast (I+{𝒘}_2{𝒆}_2^t)\dots (I+{𝒘}_{k-1}{𝒆}_{k-1}^t)$

$=I+\sum _{i=1}^{n-1}{𝒘}_i{𝒆}_i^t$

L'ultimo passaggio vale perché, ad esempio

$(I+{𝒘}_1{𝒆}_1^t)(I+{𝒘}_2{𝒆}_2^t)=I+{𝒘}_1{𝒆}_1^t+{𝒘}_2{𝒆}_2^t+{𝒘}_1{𝒆}_1^t{𝒘}_2{𝒆}_2^t$

${𝒆}_1^t{𝒘}_2=0$ come prima, perché hanno zeri in posizioni complementari. In generale

${𝒆}_i^t{𝒘}_k=0\mathrm{\forall }i\le k$

i termini misti si annullano, e quindi rimane

$(I+{𝒘}_1{𝒆}_1^t)(I+{𝒘}_2{𝒆}_2^t)=I+{𝒘}_1{𝒆}_1^t+{𝒘}_2{𝒆}_2^t$

e, procedendo così fino al passo $n$-esimo, le colonne ${𝒘}_k{𝒆}_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𝒚=𝒃$. 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𝒙=𝒃$, 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𝒙=𝒃$ al sistema $U𝒙=𝒄$ con $𝒄$ soluzione di $L𝒚=𝒃$.

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 ${ℂ}_{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 {ℂ}_{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$.