4
Fattorizzazione LU
È una fattorizzazione di una matrice in una matrice triangolare inferiore L, una
matrice triangolare superiore U. 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.
L’idea di fondo è che il metodo di eliminazione di Gauss, che permette di invertire
una matrice, può essere visto come un caso speciale della fattorizzazione LU.
Tuttavia, l’inversione diretta di una matrice è generalmente poco efficiente e
soggetta a errori numerici, mentre la fattorizzazione LU offre un approccio più
stabile e computazionalmente vantaggioso, soprattutto quando bisogna risolvere
più sistemi con la stessa matrice A.
Esempio
1. scrivo la matrice A
2 3
A=( )
4 7
2. Vado a definire L e U
per farlo mi creo una matrice triangolare inferiore L
i−1
1
lij =
(aij − ∑ lik ukj ), per i > j
ujj
k=1
1 0
L=( )
2 1
e la matrice triangolare superiore U
4 1
i−1
uij = aij − ∑ lik ukj ,
per j ≥
k=1
2 3
U =( )
0 1
3. Vado a fare LU = A
vado a fare il prodotto riga x colonna tra L e U in modo tale da ottenre A
1 0 2 3 2 3
LU = A => ( )( )=( )
2 1 0 1 4 7
I calcoli sono:
4 2
Formalmente, sia A una matrice quadrata invertibile. Essa può essere decomposta
nel prodotto di due matrici più semplici: una matrice triangolare inferiore L (con 1
sulla diagonale principale) e una matrice triangolare superiore U. Questa scrittura
prende la forma:
A = LU
e prende il nome di fattorizzazione LU.
L’idea alla base di questa decomposizione è molto simile al procedimento
dell’eliminazione di Gauss. Infatti, nell’eliminazione gaussiana si cerca di risolvere
un sistema lineare del tipo:
Ax = b
trasformando la matrice A in una forma triangolare superiore U. In questo modo il
sistema diventa Ux=b′, che può essere risolto facilmente tramite sostituzione
all’indietro.
La differenza principale è che, nella decomposizione LU, non si applica la
trasformazione anche al vettore bbb: si lavora solo sulla matrice A. Possiamo
quindi scrivere:
Ax = LUx = b
Questo approccio è molto utile, perché una volta trovata la fattorizzazione LU di A,
la si può riutilizzare per risolvere più sistemi lineari con la stessa matrice A ma con
vettori b diversi, risparmiando tempo di calcolo.
Per risolvere concretamente il sistema con la fattorizzazione LU, si introduce un
vettore intermedio y e si procede in due fasi:
Ly = b (si risolve con sostituzione in avanti)
Ux = y (si risolve con sostituzione all’indietro)
Risolvendoli in sequenza, si ottiene il vettore x cercato.
Fattorizzazione LU e determinante
4 3
La fattorizzazione A=LU è molto utile per calcolare il determinante di una matrice
in modo efficiente.
Infatti, il determinante di un prodotto di matrici è il prodotto dei determinanti:
det(A) = det(L) ⋅ det(U )
Nella fattorizzazione LU, la matrice L è triangolare inferiore con 1 sulla diagonale,
quindi
det(L) = 1
La matrice U è triangolare superiore, quindi il suo determinante è il prodotto degli
elementi sulla diagonale principale:
det(U ) = u11 ⋅ u22 ⋯ unn
Quindi, in modo compatto:
n
A = LU => det(A) = det(L) ⋅ det(U ) = 1 ⋅ ∏ = Ui
i=1
Collegamento invertibilita’
Una matrice è invertibile se e solo se il suo determinante è diverso da zero.
Durante la fattorizzazione LU, questo significa che nessun elemento diagonale di
U deve essere zero (altrimenti il determinante si annullerebbe).
In pratica:
se det(A) =
0→ A è invertibile → la fattorizzazione LU è possibile;
se det(A) = 0→ A non è invertibile → la fattorizzazione LU non esiste (o
fallisce senza pivoting).
Inoltre, una volta trovata la decomposizione A=LU, si può calcolare A−1
risolvendo n sistemi del tipo Axi = ei , cioè una colonna alla volta, usando
forward e back substitution
4 4
Sistemi multipli con LU
Quando si scrive una sequenza di equazioni del tipo
y 1 = b 1 , U x1 = y 1
Ly2 = b2 , U x2 = y2
fino a
Lym = bm , U xm = ym
si sta rappresentando una serie di sistemi lineari che hanno la stessa matrice dei
coefficienti A, ma termini noti diversi b1 , b2 , … , bm .
In pratica, invece di risolvere ogni sistema Axi = bi da zero, si effettua una sola
volta la fattorizzazione della matrice A come A=LU, dove L è triangolare inferiore e
U è triangolare superiore. Una volta ottenute L e U, per ogni vettore bi si procede
in due passaggi successivi:
prima si risolve il sistema triangolare inferiore Lyi = bi con la sostituzione in
avanti, ottenendo il vettore intermedio yi
poi si risolve il sistema triangolare superiore U xi = yi con la sostituzione
all’indietro, trovando la soluzione finale xi .
Questa scrittura con indici da 1 a m serve dunque a indicare che il procedimento si
ripete più volte, una per ogni vettore dei termini noti, ma usando sempre la stessa
decomposizione LU. Ciò è molto vantaggioso dal punto di vista computazionale,
perché la fattorizzazione viene calcolata una sola volta e poi riutilizzata per
risolvere rapidamente tutti i sistemi.
Questo approccio è particolarmente utile, ad esempio, quando si deve calcolare
l’inversa di una matrice: in quel caso, ciascun bib_ibi corrisponde a una colonna
della matrice identità, e le soluzioni xi formano le colonne dell’inversa. In
generale, la notazione Lyi = bi , U xi = yi rappresenta quindi un modo compatto
per indicare la risoluzione efficiente di più sistemi lineari che condividono la stessa
matrice ma hanno dati diversi.
Fattorizzazione LU con pivoting
4 5
Nel caso ideale, una matrice A può essere decomposta come A=LU, dove L è una
matrice triangolare inferiore e U è una matrice triangolare superiore. Tuttavia, nella
pratica numerica questo non sempre è possibile in modo stabile: durante
l’eliminazione di Gauss, può accadere che alcuni pivot (cioè gli elementi sulla
diagonale principale usati per l’eliminazione) siano molto piccoli o addirittura nulli.
Per evitare instabilità e divisioni per numeri troppo piccoli, si introduce una tecnica
chiamata pivoting, che consiste nello scambiare tra loro alcune righe della matrice
per scegliere pivot più grandi e migliorare la precisione del calcolo.
Quando si applica il pivoting, l’ordine delle righe di A viene modificato; di
conseguenza, la decomposizione non è più semplicemente A=LU, ma diventa:
P A = LU
dove P è una matrice di permutazione, cioè una matrice quadrata che rappresenta
gli scambi di righe effettuati. In altre parole, P moltiplicata per A produce una
matrice A’ in cui le righe sono state riordinate per garantire la stabilità
dell’eliminazione:
A′ = P A
La fattorizzazione corretta è quindi PA=LU, e non più A=LU.
Quando dobbiamo risolvere un sistema lineare Ax=b, l’introduzione della matrice
di permutazione modifica il procedimento. Infatti, moltiplicando entrambi i membri
per P, otteniamo:
P Ax = P b
Ponendo b′ = P b, possiamo riscrivere il sistema come:
LUx = b′
A questo punto, si procede come nel caso standard: si risolve prima Ly = b′
tramite sostituzione in avanti, e poi Ux = y tramite sostituzione all’indietro.
In sintesi, l’aggiunta del pivoting non cambia la logica generale della
fattorizzazione LU, ma introduce la matrice di permutazione P, necessaria per
4 6
mantenere la stabilità numerica e garantire che la decomposizione rimanga
corretta. In questo modo, la relazione finale è:
P A = LU
che equivale a dire che LU=A, a patto di considerare anche gli scambi di righe
rappresentati da P.
Esempio
Consideriamo:
0 2 1
A= 1 1 0
2 1 3
Notiamo subito che il primo elemento a11 = 0non può essere usato come pivot,
perché non si può dividere per zero. Qui entra in gioco il pivoting, quindi dobbiamo
scambiare la prima riga con una riga con un elemento maggiore nella prima
colonna.
Il massimo nella prima colonna è 2 (riga 3), quindi scambiamo la riga 1 con la riga
3.
La matrice di permutazione che rappresenta questo scambio è:
0 0 1
P = 0 1 0
1 0 0
Moltiplicando P⋅A, otteniamo:
2 1 3
PA = 1 1 0
= A′
0 2 1
Ora possiamo procedere con la fattorizzazione LU su A′.
4 7
Primo pivot u11 = 2, calcoliamo gli elementi di Lsotto la diagonale:
1 0
l21 = = 0.5, l31 = =0
2 2
Aggiorniamo U:
2 1 3
U= 0 ? ?
0 ? ?
(continuiamo normalmente il procedimento standard di eliminazione gaussiana su
A')
Se abbiamo un vettore b, il sistema diventa:
P Ax = LUx = P b = b′
Poi risolviamo:
1. Ly=b’ con forward substitution
2. Ux=y con back substitution
Matrice di permutazione P in informatica
Dal punto di vista teorico, nella fattorizzazione con pivoting si introduce la matrice
di permutazione P per rappresentare gli scambi di righe eseguiti durante
l’eliminazione di Gauss. Questa matrice è utile per scrivere la relazione PA=LU in
forma compatta e rigorosa.
Tuttavia, dal punto di vista computazionale, non è sempre necessario costruire
esplicitamente la matrice P. In pratica, P serve solo a memorizzare l’ordine delle
righe dopo gli scambi: quindi ciò che realmente conta è conoscere quali righe
sono state scambiate, non la matrice completa.
Creare e salvare la matrice P richiederebbe spazio in memoria proporzionale a n2
(dove n è la dimensione della matrice), il che può risultare inefficiente soprattutto
per matrici di grandi dimensioni. Per questo motivo, nei programmi e nelle librerie
4 8
di algebra lineare, non si costruisce mai esplicitamente [Link] suo posto si utilizza un
vettore di indici, chiamati ipiv (da index pivot).
Il vettore ipiv memorizza semplicemente quale riga è stata scambiata con quale, in
ciascun passo dell’eliminazione. In questo modo, l’informazione contenuta in P è
preservata, ma con un costo di memoria molto inferiore.
Quindi la matrice P è concettualmente importante per la teoria, ma nella pratica
informatica viene sostituita da un vettore di permutazioni come ipiv, che svolge la
stessa funzione in modo più efficiente.
Esempio
6 3 −2
′
b= 0 ipiv = 1
b = 6
−2 2 0
Fattorizzazione LU inplace
La fattorizzazione LU in-place è una modalità di calcolo in cui la decomposizione
A=LU viene eseguita direttamente sulla matrice originale A, cioè senza creare
nuove matrici separate per L e U.
In altre parole, invece di copiare i dati di A in due nuove matrici (una per L e una
per U), gli elementi di L e U vengono memorizzati all’interno della stessa matrice
A durante il processo di eliminazione di Gauss.
Concretamente:
Gli elementi sopra la diagonale (inclusa la diagonale stessa) della matrice
risultante contengono i valori della matrice triangolare superiore U;
Gli elementi sotto la diagonale contengono i moltiplicatori usati
nell’eliminazione, cioè i coefficienti della matrice triangolare inferiore L, che
ha diagonale unitaria.
Questo approccio è molto utilizzato nei programmi di calcolo numerico perché:
riduce il consumo di memoria, evitando di creare copie della matrice;
4 9
migliora le prestazioni, poiché lavora direttamente sui dati già allocati;
mantiene comunque tutta l’informazione necessaria per risolvere il sistema
Ax=b.
Quando serve effettivamente ricostruire L e U separatamente, è possibile estrarle
dalla matrice modificata:
L si ottiene prendendo la parte sotto la diagonale e aggiungendo una
diagonale di 1;
U si ottiene prendendo la parte superiore (compresa la diagonale).
Esempio
Prendiamo la matrice
2 4 6
A= 1 3 2
3 8 1
Vogliamo fattorizzarla come A=LU, ma senza creare nuove matrici
1. Eliminazione gaussiana
Primo pivot: a11 = 2
Calcoliamo i moltiplicatori per annullare gli elementi sotto la prima riga:
1 3
m21 = = 0.5, m31 = = 1.5
2 2
Sottraiamo m21 volte la riga 1 dalla riga 2, e m31 volte la riga 1 dalla riga 3:
R2 ← R2 − 0.5 R1
R3 ← R3 − 1.5 R1
Otteniamo:
2 4 6
0 1 −1
0 2 −8
4 10
Secondo pivot: a22 = 1
Calcoliamo il moltiplicatore per la riga 3:
2
m32 = =2
1
Sottraiamo 2×riga 2 da riga 3:
R3 ← R3 − 2R2
Risultato finale:
2 4 6
A= 0 1 −1
0 0 −6
2. Interpretazione In-place:
Ora la matrice A contiene sia L che U:
2 4 6
A= 0.5 1 −1
1.5 2 −6
(Abbiamo riscritto sotto la diagonale i moltiplicatori trovati in precedenza).
Quindi:
1 0 0 2 4 6
L= 0.5 1 0
, U= 0 1 −1
1.5 2 1 0 0 −6
Condizionamento di un problema
matematico
Problema P, risoluzione di:
4 11
{
4.19x + 7y = 15 100
= ( )
2.1x + 3.5y = 8 −57, ...
Problema P’, risoluzione di:
{
4.192x + 7y = 15 125
= ( )
2.1x + 3.5y = 8 −72, ...
Una perturbazione sulla quarta cifra significativa del primo coefficiente, ha fatto
in modo che la soluzione fosse completamente sbagliata, e quindi le soluzioni
erano lontane anche se i problemi sembravano approssimativamente vicini.
Se problemi vicini hanno soluzioni vicine, si dice che il problema e’ ben
condizionato, se invece problemi vicini hanno soluzioni lontante, si diche che il
problema e’ mal condizionato
P: d=dati, s=soluzioni
P̂: d̂ = dat perturbati, ŝ=soluzione perturbata
In generale:
4 12
σ = μδ
dove
σ = errore soluzione
δ = errore dati
dipendono entrambi da μ, che e’ l’indice di condizionamento.
💡 Se round-off e’ 10−16 e μ
problema
e’ 10−20 , e’ inutile andare a risolvere il
μdipende dal modulo della derivata, e misura la sensibilità della soluzione rispetto
ai dati:
Per una funzione f(x), è legato al modulo della derivata ∣f′(x)∣.
Se ∣f′(x)∣ è grande → problema mal condizionato.
Se ∣f′(x)∣ è piccolo → problema ben condizionato.
Nel caso f (x) = x, vicino a x=0 la derivata tende a infinito →problema mal
condizionato, perché piccole variazioni in x generano grandi variazioni in f(x).
Indice di condizionamento assoluto e
relativo
Assoluto:
CA (f , x) = ∣f ′ (x)∣
Relativo
f ′ (x)
CR (f , x) =
f (x)
4 13
Collegamento tra determinante e μ
Il det =
0⇒ e’ invertibile ed ha rango massimo, ed e’ una cosa buona, perche’
1
μ=
det =
0
ed avendo rango massimo, tutte le righe della matrice indipendenti tra loro e
quindi una relazione tra i dati migliore
Invece se det = 0
1
μ= ≈∞
0
e’ il contrario
Numero condizionamento di una matrice
Il numero di condizionamento di una matrice serve a capire quanto la soluzione di
un sistema lineare è sensibile agli errori o alle piccole variazioni nei dati di
partenza. In pratica, indica quanto un piccolo errore nel vettore dei dati bbb o nei
coefficienti della matrice A può amplificarsi nella soluzione x.
Quando si risolve un sistema Ax=b, la soluzione si ottiene formalmente come x =
A−1 b. Se nel calcolo o nei dati compare una piccola perturbazione, questa viene
amplificata dall’inversa di A. Per misurare quanto questa amplificazione può
essere grande, si introducono le norme, che servono per quantificare la
dimensione o la lunghezza di vettori e matrici.
La norma della matrice A, indicata con || A ||, rappresenta quanto la matrice può
allungare un vettore in ingresso. Allo stesso modo, la norma della matrice inversa
∣∣A−1 ∣∣misura quanto un errore può essere amplificato quando si cerca di risalire
alla soluzione. Il prodotto di queste due quantità,
μ(A) = ∣∣A∣∣ ⋅ ∣∣A−1 ∣∣
si chiama numero di condizionamento della matrice.
4 14
Questo valore è sempre maggiore o uguale a 1. Se il numero di condizionamento è
vicino a 1, il problema si dice ben condizionato, perché piccole variazioni nei dati
producono solo piccoli cambiamenti nella soluzione. Se invece il numero di
condizionamento è molto grande, il problema è mal condizionato, perché anche
un errore minimo nei dati può causare un errore molto grande nella soluzione.
4 15