Post on 16-Jul-2020
Fattorizzazione LU ed eliminazione gaussiana
Alvise Sommariva
Universita degli Studi di PadovaDipartimento di Matematica
3 maggio 2015
Alvise Sommariva Fattorizzazione LU ed eliminazione gaussiana 1/ 29
Introduzione
Problema. (Sistema lineare)
Sia A ∈ Rn×n una matrice a coeff. reali, b ∈ Rn un vettorecolonna e supponiamo di dover calcolare un vettore colonnax∗ ∈ Rn cosicche
A · x∗ = b.
Nota.
Come e noto questo problema ha soluzione unica x∗ se e solo se
det (A) 6= 0
(matrice non singolare). Ci porremo di seguito in queste ipotesi.
Alvise Sommariva Fattorizzazione LU ed eliminazione gaussiana 2/ 29
Matrici di permutazione
Definizione (Matrice di permutazione)
Una matrice P si dice di permutazione se si ottiene permutandole righe della matrice identica I .
Esempio
In questo esempio, P e ottenuta da I scambiando la prima riga conla seconda.
I =
1 0 00 1 00 0 1
, P =
0 1 01 0 00 0 1
Alvise Sommariva Fattorizzazione LU ed eliminazione gaussiana 3/ 29
Proprieta della matrice di permutazione
Nota.
Se P ∈ Rn×n si ottiene dalla matrice identica In ∈ Rn×n
scambiando la i-sima riga con la j(i)-sima allora la matriceB = PA si ottiene da A scambiando la j(i)-sima riga con lai-sima.
La matrice P ∈ Rn×n e unitaria cioe PPT = PTP = In.
Essendo PPT = PTP = In, la matrice P e invertibile e hainversa P−1 = PT . Cio implica che det(P) 6= 0.
1 Se A = BC per il teorema di Binet, alloradet(A) = det(B) det(C );
2 det(In) = 1;3 det(P) = det(PT ).
Quindi
1 = det(In) = det(PPT ) = det(P) det(PT ) = (det(P))2 ⇒ det(P) = ±1.Alvise Sommariva Fattorizzazione LU ed eliminazione gaussiana 4/ 29
Effetto della matrice di permutazione
>> % E f f e t t o d i P su A i n PA .>> A=[1 2 3 ; 4 5 6 ; 7 8 9 ]A =
1 2 34 5 67 8 9
>> P=[0 1 0 ; 1 0 0 ; 0 0 1 ]P =
0 1 01 0 00 0 1
>> P∗Aans =
4 5 61 2 37 8 9
>> %Vediamo a d e s s o che P e ’ u n i t a r i a .>> P∗P ’ans =
1 0 00 1 00 0 1
>> P ’∗ P
ans =1 0 00 1 00 0 1
>> d e t ( P )ans =−1
>>Alvise Sommariva Fattorizzazione LU ed eliminazione gaussiana 5/ 29
Matrici triangolari
Definizione (Matrici triangolari)
Una matrice A = (ai ,j) si dice
triangolare superiore, se ai ,j = 0 per i > j ;
triangolare inferiore, se ai ,j = 0 per i < j .
>> A=[1 2 3 ; 4 5 6 ; 7 8 9 ]A =
1 2 34 5 67 8 9
>> L= t r i l ( A ) % TRIANGOLARE INFERIOREL =
1 0 04 5 07 8 9
>> U=t r i u ( A ) % TRIANGOLARE SUPERIOREU =
1 2 30 5 60 0 9
>>
Alvise Sommariva Fattorizzazione LU ed eliminazione gaussiana 6/ 29
Fattorizzazione LU
Problema. (Fattorizzazione LU)
Sia A ∈ Rn×n. Determinare, se esistono,
L = li ,j triangolare inferiore con elementi diagonali uguali a 1,cioe li ,i = 1,
U triangolare superiore,
cosicche
A = LU.
Nota.
Si osservi che si chiede A = LU e non A = L + U.
Alvise Sommariva Fattorizzazione LU ed eliminazione gaussiana 7/ 29
Fattorizzazione LU
Teorema (Fattorizzazione LU e submatrici principali)
Sia A ∈ Rn×n.
Si supponga che tutte le sottomatrici principali di testaA(k) = (ai ,j)i ,j=1,...,k , k = 1, . . . , n − 1 siano non singolari.
Allora esiste ed e unica la fattorizzazione LU di A.
Nota. (Controesempio)
Non tutte le matrici posseggono la fattorizzazione LU. Un esempioin cui non esistono tali L, U per cui A = LU e la matrice
A =
(0 11 0
).
Alvise Sommariva Fattorizzazione LU ed eliminazione gaussiana 8/ 29
Fattorizzazione PA=LU
Teorema (Fattorizzazione PA=LU)
Sia A ∈ Rn×n. Allora esiste una matrice di permutazione P taleche PA = LU.Di conseguenza
La fattorizzazione A = LU non e sempre possibile.La fattorizzazione PA = LU e sempre possibile.
Nota. (Pivoting)
Per determinare la fattorizzazione PA = LU si usa una variantedell’algoritmo che determina A = LU (se esistente), ma che utilizzala tecnica del pivoting. Per dettagli, si veda [1, p.511], [2, p.172].
Nota. (Storia)
Il metodo di Gauss e stato uno dei primi implementati su uncalcolatore (1947). Obiettivo: risolvere un sistema lineare 8 per 8.
Alvise Sommariva Fattorizzazione LU ed eliminazione gaussiana 9/ 29
Fattorizzazione PA=LU in Matlab
Vediamo di seguito come eseguire la fattorizzazione PA = LU inMatlab.>> h e l p l u
l u l u factorization .[ L , U ] = l u ( A ) stores an upper triangular matrix in U and a
”psychologically l o w e r triangular matrix” ( i . e . a product of l o w e rtriangular and permutation matrices ) in L , so that A = L∗U . A can be
rectangular .
[ L , U , P ] = l u ( A ) returns unit l o w e r triangular matrix L , uppertriangular matrix U , and permutation matrix P so that P∗A = L∗U .
. . .
>>
Dall’help si capisce che
[L, U, P] = lu(A)
produce la fattorizzazione desiderata.
Alvise Sommariva Fattorizzazione LU ed eliminazione gaussiana 10/ 29
Fattorizzazione PA=LU in Matlab, esempio
>> A=[1 2 3 ; 4 5 6 ; 7 8 9 ]A =
1 2 34 5 67 8 9
>> % LA MATRICE E ’ SINGOLARE , RIGHE PROPORZIONALI ! !>> d e t ( A )ans =
6.6613 e−16>> [ L , U , P ]= l u ( A )L =
1.0000 0 00 .1429 1 .0000 00 .5714 0 .5000 1 .0000
U =7.0000 8 .0000 9 .0000
0 0 .8571 1 .71430 0 0 .0000
P =0 0 11 0 00 1 0
>> norm ( P∗A−L∗U )ans =
0>> % QUESTO CI DICE CHE P∗A−LU=0.
Alvise Sommariva Fattorizzazione LU ed eliminazione gaussiana 11/ 29
Fattorizzazione di Cholesky A = LLT per matricisimmetriche definite positive
Teorema (Fattorizzazione di Cholesky A = LLT per matricisimmetriche definite positive)
Sia A ∈ Rn×n una matrice
simmetrica, cioe A = AT ,
definita positiva, cioe avente tutti gli n autovalori λkstrettamente positivi, cioe λk > 0, per k = 1, . . . , n.
Allora esiste ed e unica la fattorizzazione di Cholesky
A = LLT
con L = (li ,j) matrice triangolare inferiore con elementi principalili ,i > 0 per i = 1, . . . , n.
Alvise Sommariva Fattorizzazione LU ed eliminazione gaussiana 12/ 29
Fattorizzazione di Cholesky in Matlab, esempio
>> h e l p c h o lc h o l Cholesky factorization .
c h o l ( A ) uses only the diagonal and upper triangle of A .The l o w e r triangle is assumed to be the ( complex conjugate )transpose of the upper triangle . If A is positive definite , then
R = c h o l ( A ) produces an upper triangular R so that R ’∗ R = A .If A is not positive definite , an e r r o r message is printed .
L = c h o l (A , ’ l o w e r ’ ) uses only the diagonal and the l o w e r triangle
of A to produce a l o w e r triangular L so that L∗L ’ = A . If
A is not positive definite , an e r r o r message is printed . When
A is s p a r s e , this syntax of c h o l is typically faster .. . .
>>
Una galleria di matrici la si puo trovare in Matlab con gallery>> h e l p g a l l e r y
g a l l e r y Higham test matrices .[ out1 , out2 , . . . ] = g a l l e r y ( matname , param1 , param2 , . . . )
takes matname , a string that is the name of a matrix family , and
the family ’ s i n p u t parameters . See the listing below f o r available
matrix families . Most of the functions take an i n p u t argument
that specifies the order of the matrix , and unless otherwise
stated , r e t u r n a single output .. . .minij Symmetric positive definite matrix MIN (i , j ) .. . .
>>
Alvise Sommariva Fattorizzazione LU ed eliminazione gaussiana 13/ 29
Fattorizzazione di Cholesky in Matlab, esempio
>> A=g a l l e r y ( ’ m i n i j ’ , 3 )A =
1 1 11 2 21 2 3
>> % V e r i f i c a numer ica che t u t t i g l i a u t o v a l o r i sono p o s i t i v i .>> e i g ( A )ans =
0.30800 .64315 .0489
>> % c h o l produce R t r i a n g o l a r e SUPERIORE , t a l e che A=R’∗R ! !>> R=c h o l ( A )R =
1 1 10 1 10 0 1
>> norm (A−R ’∗ R )ans =
0>> % c h o l (A, ’ lower ’ ) produce L t r i a n g o l a r e INFERIORE , t a l e che A=L∗L ’ ! !>> L=c h o l (A , ’ l o w e r ’ )L =
1 0 01 1 01 1 1
>> norm (A−L∗L ’ )ans =
0>>
Alvise Sommariva Fattorizzazione LU ed eliminazione gaussiana 14/ 29
Eliminazione Gaussiana
Si supponga di dover risolvere Ax = b con
A ∈ Rn×n, det(A) 6= 0 (cioe A non singolare),
b ∈ Rn×1 = Rn.
x∗ unica soluzione del sistema lineare Ax = b, cioe Ax∗ = b.
Se PA = LU allora essendo det(P) 6= 0, abbiamo che
Ax∗ = b ⇔ PAx∗ = Pb ⇔ LUx∗ = Pb.
Posto y∗ = Ux∗, da LUx∗ = Pb abbiamo che y∗ e lasoluzione del sistema triangolare inferiore Ly∗ = Pb.
Una volta ottenuto y∗, essendo Ux∗ = y∗, x∗ e la soluzionedel sistema triangolare superiore Ux∗ = y∗.
Alvise Sommariva Fattorizzazione LU ed eliminazione gaussiana 15/ 29
Eliminazione Gaussiana
Questa osservazione suggerisce il seguente metodo per risolvereAx = b con A non singolare.
Metodo (Eliminazione gaussiana)
Si determini la fattorizzazione PA = LU di A ∈ Rn×n (costocomputazionale O(n3/3)).
Si determini la fattorizzazione c = Pb.
Si risolva il sistema triangolare inferiore Ly = c (costocomputazionale O(n2/2)).
Si risolva il sistema triangolare superiore Ux = y (costocomputazionale O(n2/2)).
Alvise Sommariva Fattorizzazione LU ed eliminazione gaussiana 16/ 29
Eliminazione Gaussiana, esempio in Matlab
>> A=g a l l e r y ( ’ m i n i j ’ , 5 )A =
1 1 1 1 11 2 2 2 21 2 3 3 31 2 3 4 41 2 3 4 5
>> x_sol=ones ( 5 , 1 ) ;>> b=A∗x_solb =
59
121415
>> % Ho un s i s t e m a Ax=b con s o l u z i o n e x s o l =[1 1 1 1 1 ] ’ .>> % APPROSSIMO LA SOLUZIONE COL METODO DI ELIMINAZIONE GAUSSIANA .>> [ L , U , P ]= l u ( A ) ; % F a t t o r i z z a z i o n e PA=LU .>> c=P∗b ; % Ax=b a l l o r a PAx=Pb . Pongo Pb=c e r i s o l v o PAx=c .>> y=L\c ; % S o l u z i o n e s i s t e m a t r i a n g o l a r e i n f e r i o r e .>> x=U\y ; % S o l u z i o n e s i s t e m a t r i a n g o l a r e s u p e r i o r e .>> x
x =11111
>> % Ho c a l c o l a t o c o r r e t t a m e n t e l a s o l u z i o n e i n quanto x s o l = x .
Alvise Sommariva Fattorizzazione LU ed eliminazione gaussiana 17/ 29
Eliminazione Gaussiana: A = LU o PA = LU?
Problema.
Se
PA = LU e ottenuta col metodo di pivoting per colonneimplementato dal Matlab,
esiste pure la fattorizzazione A = LU,
quale delle due e da preferire?
Per questioni di stabilita e da preferire il metodo tramite pivotingper colonne. Illustriamo questo in un esempio.
Alvise Sommariva Fattorizzazione LU ed eliminazione gaussiana 18/ 29
Eliminazione Gaussiana: A = LU o PA = LU?
Esempio (Matrice di Hankel)
Sia H(n) ∈ Rn×n la matrice di Hankel di ordine n, i cui elementisono definiti come segue
H(n)i ,n+k−i =
{2k se k > 0
21/(2−k) se k ≤ 0
con i = 1, . . . , n, k = i + 1− n, . . . , i . La matrice H(n) e invertibile.
Sia x∗ = [1, . . . , 1] ∈ Rn×1 e b = Ax. Ovviamente x∗ e l’unicasoluzione di Ax = b.
Sia xLU la soluzione ottenuta con il metodo di Eliminazionegaussiana senza permutazione.
Sia xLUP la soluzione ottenuta con il metodo di Eliminazionegaussiana con permutazione (dovuta a pivoting per colonne).
Alvise Sommariva Fattorizzazione LU ed eliminazione gaussiana 19/ 29
Eliminazione Gaussiana: A = LU o PA = LU?
eLU = ‖x∗ − xLU‖2 =√∑n
k=1(x∗k − xLUk )2.
eLUP = ‖x∗ − xLUP‖2 =√∑n
k=1(x∗k − xLUPk )2.
Si verifica sperimentalmente che
n eLU eLUP n eLU eLUP1 0.00e + 00 0.00e + 00 13 1.14e − 09 1.18e − 112 1.11e − 16 8.88e − 16 14 9.34e − 09 1.26e − 113 4.22e − 15 3.00e − 15 15 3.96e − 08 4.11e − 114 4.22e − 15 4.55e − 15 16 2.36e − 07 6.78e − 115 1.31e − 14 1.58e − 14 17 2.10e − 06 3.23e − 106 3.48e − 13 2.60e − 14 18 1.18e − 05 3.09e − 107 1.24e − 13 3.73e − 14 19 4.10e − 05 5.63e − 108 1.39e − 12 1.59e − 13 20 1.67e − 04 1.26e − 099 7.13e − 12 6.02e − 13 21 4.39e − 04 3.22e − 0910 1.05e − 11 3.41e − 13 22 2.16e − 02 3.80e − 0911 1.56e − 11 9.58e − 13 23 2.97e − 02 1.03e − 0812 1.57e − 10 4.13e − 12 24 3.65e − 02 1.58e − 08
Alvise Sommariva Fattorizzazione LU ed eliminazione gaussiana 20/ 29
Eliminazione Gaussiana con A simmetrica definita positiva
Si supponga di dover risolvere Ax = b con
A ∈ Rn×n, simmetrica e definita positiva.
b ∈ Rn×1 = Rn.
x∗ unica soluzione del sistema lineare Ax = b, cioe Ax∗ = b.
Se A = LLT allora abbiamo che
Ax∗ = b ⇔ LLT x∗ = b.
Posto y∗ = LT x∗, da L ∗ y∗ = LLT x∗ = b abbiamo che y∗ ela soluzione del sistema triangolare inferiore Ly∗ = b.
Una volta ottenuto y∗, essendo LT x∗ = y∗, x∗ e la soluzionedel sistema triangolare superiore LT x∗ = y∗.
Alvise Sommariva Fattorizzazione LU ed eliminazione gaussiana 21/ 29
Eliminazione Gaussiana con A simmetrica definita positiva
Questa osservazione suggerisce il seguente metodo per risolvereAx = b con A simmetrica e definita positiva.
Metodo (Eliminazione gaussiana se A simmetrica e definitapositiva)
Si determini la fattorizzazione A = LLT di A ∈ Rn×n (costocomputazionale O(n3/6)).
Si risolva il sistema triangolare inferiore Ly = c (costocomputazionale O(n2/2)).
Si risolva il sistema triangolare superiore LT x = y (costocomputazionale O(n2/2)).
Alvise Sommariva Fattorizzazione LU ed eliminazione gaussiana 22/ 29
Eliminazione Gaussiana, esempio in Matlab
>> A=g a l l e r y ( ’ m i n i j ’ , 5 )A =
1 1 1 1 11 2 2 2 21 2 3 3 31 2 3 4 41 2 3 4 5
>> ( e i g ( A ) ) ’ans =
0.2716 0 .3533 0 .5830 1 .4487 12.3435>> % A e ’ s i m m e t r i c a d e f i n i t a p o s i t i v a .>> x_sol=ones ( 5 , 1 ) ; b=A∗x_sol ;>> % Ho un s i s t e m a Ax=b con s o l u z i o n e x s o l =[1 1 1 1 1 ] ’ .>> % APPROSSIMO LA SOLUZIONE COL METODO DI ELIMINAZIONE GAUSSIANA VIA CHOLESKY .>> L=c h o l (A , ’ l o w e r ’ )L =
1 0 0 0 01 1 0 0 01 1 1 0 01 1 1 1 01 1 1 1 1
>> y=L\b ;>> x=L ’\ yx =
11111
>> % Ho c a l c o l a t o c o r r e t t a m e n t e l a s o l u z i o n e i n quanto x s o l = x .
Alvise Sommariva Fattorizzazione LU ed eliminazione gaussiana 23/ 29
Esercizi
Esercizio
Implementare una routine Matlab
flag = issymm(A).
che calcolando gli autovalori di una matrice A, stabilisca se Ae simmetrica definita positiva (se flag=1 allora e simmetricadefinita positiva altrimenti non lo e) .
Nota: per vedere che e simmetrica basta notare che cio e verose
norm(A− A′) = 0.
Alvise Sommariva Fattorizzazione LU ed eliminazione gaussiana 24/ 29
Esercizi
Esercizio
Implementare una routine Matlab
x = linear solver(A, b).
che
calcoli la soluzione x∗ mediante l’eliminazione gaussiana viafattorizzazione di Cholesky se A e simmetrica e definitapositiva,altrimenti
si effettui PA = LU ove U = (ui,j) e osservato che
det(A) = det(P) det(U) = ±n∏
k=1
uk,k
si verifichi se A e o meno singolare;se A e non singolare, si determini la soluzione x∗ conl’eliminazione gaussiana via fattorizzazione PA = LU.
Alvise Sommariva Fattorizzazione LU ed eliminazione gaussiana 25/ 29
Esercizi
Esercizio
Testare il codice precedente per risolvere il problema Ax = bdove
A = gallery(′poisson′, 20); b = ones(size(A, 1));
Testare che soluzione x ottenuta coincida con quella fornitadal Matlab via
x sol = A\b;
A tal proposito testare che
norm(x− x sol);
sia molto piccola (ad esempio dell’ordine di 10−15).
Alvise Sommariva Fattorizzazione LU ed eliminazione gaussiana 26/ 29
Esercizi
Esercizio (Facoltativo)
Si scarichino i files solve linear LUP.m, solve linear LUP.m,hankel matrix.m, dalla directory del corso.
La chiamata
x = solve linear LU(A, b);
risolve il sistema Ax = b mediante Eliminazione Gaussiana viafattorizzazione A = LU (se esiste!).
La chiamata
x = solve linear LUP(A, b);
risolve il sistema Ax = b mediante Eliminazione Gaussiana viafattorizzazione PA = LU.
Alvise Sommariva Fattorizzazione LU ed eliminazione gaussiana 27/ 29
Esercizi
La chiamataA = hankel matrix(n)
produce una matrice di Hankel A ∈ Rn×n.
Sia b = Ax dove
A=hankel matrix(20);
x=ones(20,1); b=A*x;
Si risolva il problema Ax = b mediante solve linear LU esia x1 la soluzione ottenuta.
Si risolva il problema Ax = b mediante solve linear LUP esia x2 la soluzione ottenuta.
Si calcolino gli errori norm(x-x1) e norm(x-x2). I risultatisono analoghi? Quale metodo e da preferire?
Alvise Sommariva Fattorizzazione LU ed eliminazione gaussiana 28/ 29
Bibliografia
K.E. Atkinson An introduction to Numerical Analysis, Wiley, (1989).
D. Bini, M. Capovani, O. Menchi Metodi numerici per l’algebra lineare, Zanichelli, (1988).
Alvise Sommariva Fattorizzazione LU ed eliminazione gaussiana 29/ 29