Algebra lineare

1
Algebra lineare

• G. Puppo

2
Algebra lineare

• Numero di condizionamento
• Sistemi triangolari
• Fattorizzazione LU
• Soluzione di sistemi lineari in Matlab
• Effetto del numero di condizionamento

3
Numero di condizionamento

• Per calcolare il numero di condizionamento,
  Matlab usa listruzione cond(a,p), dove
• p 1 è per la norma 1
• p 2 è per la norma 2 (default)
• p inf è per la norma infinito

4
Andamento del numero di condizionamento
Costruiamo un programma che calcoli il numero di
condizionamento al variare di N per matrici N per
N assegnate. In particolare il programma deve

• avere un ciclo su N
• calcolare le matrici richieste per ogni N
• calcolare il numero di condizionamento
  corrispondente per ogni matrice, in ognuna delle
  norme 1, 2 e inf
• stampare un grafico con i risultati

5
function plot_con(nmax) Calcola il numero di
condizionamento per le matrici RAND(N), con
N1NMAX, confrontando su un grafico i risultati
ottenuti in norma 1, norma 2, e norma inf for
n1nmax arand(n) con1(n) cond(a,1)
con2(n) cond(a,2) coninf(n)
cond(a,inf) end plot(log10(con2)) hold
on plot(log10(con1),'g') plot(log10(coninf),'m')(l
og10(
6
(No Transcript)
7
Proviamo altre matrici ...

• Cambiamo le matrici calcolate nella function
  plot_con
• Sostituiamo a rand(n) le matrici tridiagonali,
  fornite dalla function tridiag(n)
• Sostituiamo le matrici di Hilbert, hilb(n), che
  hanno un condizionamento exp(n).

8
Matrice tridiagonale tridiag(n)
9
Matrici di Hilbert
La scala verticale è logaritmica
10
Sistemi triangolari inferiori
Scrivo una function che calcola la soluzione del
sistema Axb, nel caso in cui A è triangolare
inferiore. In una function che deve essere
utilizzata da diversi utenti è bene controllare
che i dati forniti siano coerenti. In questo caso
controllo che - A sia quadrata - A sia
triangolare inferiore - A sia non singolare - Il
vettore b abbia le dimensioni giuste
11
Controllo che la matrice A sia ben impostata
function xtri_inf(a,b) TRI_INF(A,B) Risolve il
sistema triangolare inferiore AXB. Se A
non e' quadrata, o A non e' triangolare
inferiore, stampa un messaggio di
errore Controlla le dimensioni di A n,m
size(a) if m n display('A non e''
quadrata') return end for i1n for
j(i1)n if a(i,j) 0
display('A non e'' triangolare inf')
return end end end
12
Termino il controllo dei dati
Controlla le dimensioni di B if length(b) n
display('B non e'' compatibile')
return end Controlla se A e' singolare for
i1n if abs( a(i,i)) lt epsnorm(b)
display('A e'' quasi singolare') return
end end
13
Finalmente, risolvo il sistema
La formula ricorsiva è
Risolve il sistema x(1) b(1)/a(1,1) for
i2n sumb(i) for j1(i-1) sum
sum - a(i,j)x(j) end x(i)
sum/a(i,i) end
Lalgoritmo corrispondente è
14
Risolvo un sistema triangolare superiore. Controll
o che A sia ben impostata
function xtri_sup(a,b) TRI_SUP(A,B) Risolve il
sistema triangolare superiore AXB. Se A
non e' quadrata, o A non e' triangolare
superiore, stampa un messaggio di
errore Controlla le dimensioni di A n,m
size(a) if m n display('A non e''
quadrata') return end for i1n for
j1i-1 if a(i,j) 0
display('A non e'' triangolare sup')
return end end end
15
Controllo che B sia compatibile e che A non sia
singolare
Controlla le dimensioni di B if length(b) n
display('B non e'' compatibile')
return end Controlla se A e' singolare for
i1n if abs( a(i,i)) lt epsnorm(b)
display('A e'' quasi singolare') return
end end
16
Risolvo il sistema
La formula ricorsiva è
Risolve il sistema x(n) b(n)/a(n,n) for
in-1-11 sumb(i) for ji1n
sum sum - a(i,j)x(j) end x(i)
sum/a(i,i) end
E lalgoritmo corrispondente è
17
Calcola la fattorizzazione LU di una matrice
function l,u elleu(a) ELLEU(A) Calcola la
fattorizzazione LU di A, senza pivoting Sintassi
L,UELLEU(A) esce con un messaggio di
errore, se trova un pivot lt
EPSNORM(A) nornorm(a) Controlla le
dimensioni di A n,m size(a) if m n
display('A non e'' quadrata') return end
18
Calcola la fattorizzazione LU, riscrivendo su A
l-eye(n) for k1n-1 if abs(a(k,k)) lt
epsnor display('Pivot troppo piccolo')
return end m(k1n,k)
-a(k1n,k)/a(k,k) a(k1n,k1n)
a(k1n,k1n)m(k1n,k)a(k,k1n) end
19
infine, immagazzina i risultati nelle matrici L
e U
Immagazzina i risultati nelle matrici l e u l
-m u zeros(n) for i1n for jin
u(i,j) a(i,j) end end
20
Risoluzione di un sistema lineare
Con le tre functions appena costruite, posso
risolvere un sistema lineare. Devo impostare una
matrice quadrata A e un vettore B di termini
noti. Per risolvere il sistema AxB, devo

• Calcolare la fattorizzazione LU
• Risolvere il sistema triangolare inferiore Lyb
• Risolvere il sistema triangolare superiore Uxy

gtgt l,uelleu(a) gtgt ytri_inf(l,b) gtgt
xtri_sup(u,y)
21
Esempio
Risolvo il sistema Axb, per
gtgt a2 -1 0 -1 2 -1 0 -1 2 gtgt a a 2
-1 0 -1 2 -1 0 -1
2 gtgt bones(3,1)
22
Calcola la fattorizzazione LU
gtgt l,uelleu(a) l 1.0000 0
0 -0.5000 1.0000 0 0
-0.6667 1.0000 u 2.0000 -1.0000
0 0 1.5000 -1.0000 0
0 1.3333
Osservo che
gtgt lu ans 2 -1 0 -1 2
-1 0 -1 2
23
Risolvo i due sistemi triangolari
gtgt ytri_inf(l,b) y 1.0000 1.5000
2.0000 gtgt xtri_sup(u,y) x 1.5000 2.0000
1.5000
Osservo che
gtgt ax'-b ans 1.0e-015 0
0.2220 -0.4441
Quindi, il residuo è piccolo, anche se non è
zero. La soluzione è accettabile
24
Provo un altro esempio con ahilb(5),
bones(5,1). Ottengo
gtgt ahilb(5) gtgt bones(5,1) gtgt
l,uelleu(a) gtgt ytri_inf(l,b) gtgt
xtri_sup(u,y) gtgt norm(ax'-b) ans
3.1776e-014
Quindi il residuo è aumentato, infatti
gtgt cond(a) ans 4.7661e005
25
Leffetto del numero di condizionamento si vede
meglio con questo esempio
Scelgo b, in modo da conoscere la soluzione
esatta, x. Calcolo la soluzione x1 risolvendo il
sistema lineare. In aritmetica esatta, avrei
x-x10. In aritmetica floating point
gtgt norm(x-x1')/norm(x) ans 1.5531e-012
gtgt ahilb(5) gtgt xones(5,1) gtgt bax gtgt
l,uelleu(a) gtgt ytri_inf(l,b) gtgt
x1tri_sup(u,y)
La differenza x-x1 è molto più grande del
residuo
26
Risoluzione di sistemi lineari con Matlab
Matlab risolve il sistema lineare Axb con il
comando xA\b.
Se A è quadrata, questo comando implica i
seguenti passi - Calcola la fattorizzazione LU
con pivoting sulle colonne - Risolvi i due
sistemi triangolari
Se A è rettangolare (sistemi sovradeterminati),
Matlab calcola la soluzione nel senso dei minimi
quadrati, cioè - Calcola la fattorizzazione QR
con pivoting sulle colonne - Risolve il sistema
triangolare superiore
Se il condizionamento di A è elevato, stampa un
messaggio di warning, ma calcola ugualmente la
soluzione.
27
Risoluzione di sistemi lineari con Matlab 2
Nella versione 7, è stata inserita una nuova
function per risolvere sistemi lineari con
struttura particolare
xlinsolve(a,b)
28
Esempio
gtgt a2 -1 0 -1 2 -1 0 -1 2 gtgt
bones(3,1) gtgt xa\b x 1.5000 2.0000
1.5000
Attenzione! Matlab usa sia / che \ ma i due
operatori hanno un effetto diverso. Provare per
credere!
29
Attenzione! Matlab calcola una soluzione anche
quando il sistema non ammette soluzione. Nellese
mpio seguente, a è singolare
gtgt a1 2 3 4 5 6 7 8 9 gtgt b1 1 0 gtgt
xa\b Warning Matrix is close to singular or
badly scaled. Results may be inaccurate.
RCOND 1.541976e-018. gtgt norm(ax-b) ans
5.0990
RCOND è il reciproco del numero di
condizionamento se RCOND è piccolo, la matrice è
mal condizionata.
30
Fattorizzazione LU
Matlab calcola la fattorizzazione LU di una
matrice con il comando
l,u lu(a)
In questo caso, l contiene anche gli scambi di
riga effettuati. Se la fattorizzazione è PA LU,
allora lP-1 L.
Oppure posso usare tre argomenti in output
l,u,plu(a)
In questo caso, lL, uU, pP.
31
Esempio
Notare che l contiene la matrice triangolare
L con le righe scambiate, cioè lP-1 L
gtgt a1 2 3 4 5 6 7 8 9 gtgt l,ulu(a) l
0.1429 1.0000 0 0.5714 0.5000
1.0000 1.0000 0 0 u
7.0000 8.0000 9.0000 0 0.8571
1.7143 0 0 0.0000
32
Per forzare la soluzione di un sistema tramite
fattorizzazione LU, devo quindi dare i comandi
gtgt l,ulu(a) gtgt yl\b gtgt xu\y
33
Fattorizzazione QR
Matlab esegue la fattorizzazione AQR mediante il
comando
gtgtq,r qr(a)
Per forzare la soluzione di un sistema lineare
mediante fattorizzazione QR devo quindi dare i
comandi
gtgt q,rqr(a) gtgt xr\q'b
34
Esercizio
35
Effetto del numero di condizionamento

• Variazione del residuo per matrici mal
  condizionate
• Effetto del condizionamento sulla precisione
  della soluzione numerica
• Stima del rango

36
Variazione del residuo per matrici mal
condizionate
Scrivo un programma che calcoli, per n da 1 a
100 - ahilb(n) - brand(n,1) - risolve il
sistema axb - calcola la norma del residuo
ax-b - calcola il condizionamento di a - stampa
un grafico con il condizionamento ed il residuo
in funzione di n
Vorrei studiare come varia il residuo r
Ax-b, in funzione del condizionamento della
matrice A. Scelgo come matrice mal condizionata
la matrice di Hilbert, mentre per b assegno un
vettore di numeri casuali
37
Script res_hilb.m
Esegue un grafico del residuo, in funzione di n,
per sistemi lineari del tipo Hilb(n)xrand(n,1)
nmax100 for n1nmax ahilb(n)
brand(n,1) xa\b c(n)cond(a)
r(n)norm(ax-b) end plot(log10(c)) hold
on plot(log10(reps),'m')
38
Ottengo questi risultati
39
Effetto del condizionamento sulla precisione
della soluzione
Scrivo un programma che calcoli, per n da 1 a
100 - ahilb(n) - xexaones(n,1)
(xexasoluzione esatta) - baxexa (b
termine noto corrispondente a xexa) - risolve il
sistema axb (xsoluzione approssimata) -
calcola lerrore relativo fra x e xexa - calcola
il condizionamento di a - stampa un grafico con
il condizionamento e lerrore relativo
Imposto un problema Axb di cui conosco la
soluzione esatta. Risolvo numericamente il
problema, calcolando la soluzione stimata y.
Quindi calcolo lerrore x-y. Mi aspetto
errori grandi se uso per A una matrice mal
condizionata, come, di nuovo, la matrice di
Hilbert.
40
Script hilb_exa.m
Calcola la differenza fra soluzione esatta e
soluzione approssimata, per i sistemi lineari
hilb(n)xb e ne stampa il grafico,
confrontando con il condizionamento nmax100 for
n1nmax ahilb(n) xexaones(n,1)
baxexa xa\b c(n)cond(a)
err(n)norm(x-xexa)/norm(xexa) end plot(log10(c))
hold on plot(log10(erreps),'m')
41
Ottengo questi risultati
42
Stima del rango di una matrice
Per stimare il rango di una matrice posso

• Calcolare la fattorizzazione LU e contare gli
  elementi U(i,i) che hanno modulo maggiore di una
  tolleranza assegnata
• Calcolare la fattorizzazione QR e contare gli
  elementi R(i,i) che hanno modulo maggiore di una
  tolleranza assegnata

Per matrici mal condizionate il secondo metodo è
più affidabile. La tolleranza assegnata ha un
ruolo fondamentale
43
Listato della function rango.m
function rrango(a,tol) Calcola il rango della
matrice A, usando il metodo QR Sintassi
RRANGO(A) Calcola il rango di A. Usa EPS
se r(i,i) ltEPSnorm(a,1) gt r(i,i)0
RRANGO(A,TOL) Calcola il rango di A,
se r(i,i) ltTOL gt r(i,i)0
Questa function può essere chiamata in due
modi rango(a) usa una tolleranza di default
TOLEPSNORM(a,1) rango(a,tol) usa la tolleranza
TOL fissata in input
44
Questa possibilità di scelta è implementata con
la function NARGIN, che conta il numero di
argomenti in input
if nargin lt 2 tolepsnorm(a,1) end
Quindi, se NARGIN ritorna il valore 1, la
function RANGO calcola la tolleranza,
TOLEPSNORM(A,1) altrimenti TOL non viene
calcolata, ma viene usato il valore assegnato in
input
45
function rrango(a,tol) Calcola il rango della
matrice A, usando il metodo QR Sintassi
RRANGO(A) Calcola il rango di A. Usa EPS
se r(i,i) ltEPSnorm(a,1) gt r(i,i)0
RRANGO(A,TOL) Calcola il rango di A,
se r(i,i) ltTOL gt r(i,i)0 if nargin lt 2
tolepsnorm(a,1) end q,rrqr(a) r0
m,nsize(a) for i1n if abs(rr(i,i))gttol
rr1 end end
46
Provo ora a calcolare il rango delle matrici di
Hilbert. Uso le seguenti istruzioni
for i1100 r(i)rango(hilb(i)) end gtgt plot(r)
Notare che tutte le matrici di Hilbert sono non
singolari, quindi dovrei avere r(i)i
47
Rango delle matrici di Hilbert, con
TOLEPSNORM(A,1)
48
Rango delle matrici di Hilbert, con TOL1e-40
Ora tutte le matrici hanno rango pieno
49
Apparentemente, i risultati ottenuti con
TOL1e-40 sono più corretti di quelli ottenuti
con TOLEPSNORM(A,1). Calcoliamo però qual è
laccuratezza della fattorizzazione QR che stiamo
calcolando, per esempio per HILB(50)
gtgt ahilb(50) gtgt q,rqr(a) gtgt norm(qr-a) ans
5.7905e-016
Quindi non ha senso usare un valore di TOL più
basso di circa 1e-16, perché la matrice R non è
abbastanza accurata