Post on 16-Feb-2019
SISTEMI LINEARI QUADRATI
Dipartimento di Sc. Matematiche ed Informatiche, Universita di Siena
LABORATORIO di CALCOLO NUMERICO
a.a. 2007–2008
SISTEMI LINEARI QUADRATI – p.1/47
RICHIAMI di ALGEBRA LINEAREDEFINIZIONI
A ∈ Rn×n simmetrica se A = AT ;
A ∈ Cn×n Hermitiana se A = A∗, dove A∗ = (aj,i) i=1,...,nj=1,...,n
;
A ∈ Rn×n ortogonale se AT A = I;
A ∈ Cn×n (simmetrica) unitaria se A∗A = I;
A ∈ Rn×n (simmetrica) definita positiva se xT Ax > 0, ∀x 6= 0;
A ∈ Cn×n (hermitiana) definita positiva se x∗Ax > 0, ∀x 6= 0;
A diagonale dominante in senso forte[debole] se
|ai,i| >
n∑
j=1j 6=i
|ai,j |,
|ai,i| ≥
n∑
j=1j 6=i
|ai,j |
i = 1, . . . , n
A diagonale dominante se
|ai,i| ≥n∑
j=1j 6=i
|ai,j |, i = 1, . . . , n, ed ∃k : |ak,k| >
n∑
j=1j 6=k
|ak,j |
SISTEMI LINEARI QUADRATI – p.2/47
RICHIAMI di ALGEBRA LINEARETEOREMASe A ∈ Rn×n è simmetrica, sono equivalenti:
1. A definita positiva;
2. gli autovalori di A sono reali e positivi;
3. i determinanti dei minori principali di A verificano: det(Ak) > 0, k = 1, . . . , n
(criterio di Sylvester).
L’algebra lineare numerica richiede spazi lineari normati
ESEMPIO: Si consideri il sistema lineare
8 1 1
1 5 −1
1 −1 5
x1
x2
x3
=
26
7
7
la cui soluzione e x = [3, 1, 1]T . Si supponga di aver ottenuto le seguenti approssimazioni della
soluzione: v = [2.9997675, 1.0002365, 1.0003005]T
z = [3.0002725, 0.9996165, 0.9995525]T
Come possiamo stabilire quale delle due grandezze e piu vicina alla soluzione del sistema?
SISTEMI LINEARI QUADRATI – p.3/47
Norma di vettoreDefinizione: Sia X spazio lineare e f : X → R funzionale tale che
1. f(x) ≥ 0, ∀x ∈ X, f(x) = 0 ⇔ x = 0;
2. f(x + y) ≤ f(x) + f(y), ∀x, y ∈ X;
3. f(αx) = |α|f(x), ∀α ∈ R(C), ∀x ∈ X;
X è detto spazio lineare normato (sln) e f è detta norma.In particolare per X = Rn si usano le norme p o Hölderiane:
‖x‖p := (|x1|p + |x2|p + · · · + |xn|p)1p
‖x‖1 := |x1| + |x2| + · · · + |xn|
‖x‖2 :=√
x21 + x2
2 + · · · + x2n
‖x‖∞ := max{|x1|, |x2|, . . . , |xn|}
SISTEMI LINEARI QUADRATI – p.4/47
Norma di vettore
Definizione: Siano f ed g due norme di Rn. Se esistono c, C ∈ R+ t.c.
cg(x) ≤ f(x) ≤ Cg(x), ∀x ∈ Rn
f ed g si dicono norme equivalenti.
‖ · ‖1, ‖ · ‖2, ‖ · ‖∞ sono norme equivalenti
‖x‖∞ ≤ ‖x‖1 ≤ n‖x‖∞‖x‖∞ ≤ ‖x‖2 ≤
√n‖x‖∞
‖x‖1/√
n ≤ ‖x‖2 ≤ ‖x‖1
ma con diversa geometria: S = {x ∈ Rn : ‖x‖ ≤ 1}
|| ||1
|| ||2
|| ||∞
SISTEMI LINEARI QUADRATI – p.5/47
Norma di matrici in Rn×n
Definizione: Sia ‖ · ‖Rn una norma vettoriale. ‖ · ‖Rn×n si dice compatibile o
consistente con ‖ · ‖Rn se ∀A, B ∈ Rn×n e ∀x ∈ Rn si ha
‖Ax‖Rn ≤ ‖A‖Rn×n‖x‖Rn
‖AB‖Rn×n ≤ ‖A‖
Rn×n‖B‖Rn×n
Definizione: Sia ‖ · ‖Rn×n compatibile, si dice indotta da ‖ · ‖Rn se ∀A ∈ Rn×n
∃x ∈ Rn t.c.
‖Ax‖Rn = ‖A‖Rn×n‖x‖Rn
Definizione: Sia data ‖ · ‖Rn norma vettoriale, il funzionale ‖ · ‖ : Rn×n → R t.c.
‖A‖ := supx6=0
‖Ax‖Rn
‖x‖Rn
= max‖x‖Rn=1
‖Ax‖Rn
è una norma, norma naturale
Ogni norma naturale è indotta:sia y : ‖Ay‖ = max‖x‖=1 ‖Ax‖ ⇒ ‖Ay‖ = ‖A‖ = ‖A‖‖y‖
SISTEMI LINEARI QUADRATI – p.6/47
Principali norme indotte naturali
‖ ‖1 induce ‖A‖1 := max1≤j≤n
n∑
i=1
|ai,j |
‖ ‖∞ induce ‖A‖∞ := max1≤i≤n
n∑
j=1
|ai,j |
‖ ‖2 induce ‖A‖2 :=√
ρ(AT A)
dove data M ∈ Rp×p, ρ(M) = max1≤k≤p |λk(M)| è il suo raggio spettrale
SISTEMI LINEARI QUADRATI – p.7/47
Proprietà delle norme naturali
Sia M ∈ Rp×p
ρ(M) = inf{‖M‖ : ∀‖ · ‖ norma naturale}
ovvero ∀ε > 0, ∃‖ · ‖ norma naturale t.c ρ(M) ≤ ‖M‖ ≤ ρ(M) + ε, ne segue
∀M ∈ Rp×p, ∀‖ · ‖ norma naturale,
ρ(M) ≤ ‖M‖
Una matrice M ∈ Rp×p si dice convergente se limk→∞ Mk = 0
Teorema
M convergente ⇔ ρ(M) ≤ 1
Invarianza della norma rispetto a trasformazioni ortogonali:
x ∈ Rn, U ∈ Rn×n ortogonale: ‖x‖2 = ‖Ux‖2
A ∈ Rm×n, U ∈ Rm×m, V ∈ Rn×n ortogonali: ‖A‖2 = ‖UAV ‖2
SISTEMI LINEARI QUADRATI – p.8/47
SISTEMI LINEARI QUADRATI
Sia A ∈ Rn×n, b ∈ Rn, se det A 6= 0 allora ∃!x ∈ Rn t.c.
Ax = b.
Molti modelli matematici significativi sono di tipo lineare
Sistemi lineari nascono da contesti diversi, è importante disporre di un ampiavarietà di algoritmi in modo da scegliere il più adatto al problema specifico(stabilità, occupazione di memoria, velocità)
I metodi per i sistemi lineari si dividono in due gruppi:
Metodi diretti: si basano sull’idea di trasformare il sistema attraverso unnumero finito di operazioni in un sistema equivalente di cui siaesplicitamente calcolabile la soluzione. In assenza di errori diarrotondamento forniscono la soluzione esatta.
Metodi iterativi: la soluzione è ottenuta come limite di una successione.Permettono di sfruttare la sparsità della matrice dei coeffcienti in quanto alcontrario dei metodi diretti, tale matrice non viene modificata
SISTEMI LINEARI QUADRATI – p.9/47
Condizionamento
Dato Ax = b si esprimono le perturbazioni sui dati con εF, F ∈ Rn×n e εf, f ∈ Rn.Si risolve il sistema
(A + εF )x(ε) = b + εf ; (1)
con x(0) = x. Se det A 6= 0 allora per ε “piccolo” det(A + εF ) 6= 0; posso ricavare
x(ε) = (A + εF )−1(b + εf)
funzione di ε derivabile in un intorno dell’origine.Derivando (1) si ha
Fx(ε) + (A + εF )x(ε) = f ;
per ε = 0 : Fx(0) + Ax(0) = f da cui x(0) = A−1(f − Fx(0)).Sviluppando in serie x(ε) = x(0) + εx(0) + O(ε2) ⇒ x(ε) − x(0) ' εx(0) da cui
‖x(ε) − x(0)‖‖x(0)‖ ' ‖εx(0)‖
‖x(0)‖ = ε‖A−1(f − Fx(0))‖
‖x(0)‖
SISTEMI LINEARI QUADRATI – p.10/47
Condizionamento
‖x(ε) − x(0)‖‖x(0)‖ ≤ ε
‖A−1‖‖f − Fx(0)‖‖x(0)‖
≤ ε‖A−1‖( ‖f‖‖x(0)‖ + ‖F‖
)
= ε‖A−1‖‖A‖( ‖f‖‖A‖‖x‖ +
‖F‖‖A‖
)
≤ ‖A−1‖‖A‖(‖εf‖
‖b‖ +‖εF‖‖A‖
)
≤ κ(A)
(‖εf‖‖b‖ +
‖εF‖‖A‖
)
κ(A) := ‖A−1‖‖A‖: numero di condizionamento (‖ · ‖ qualsiasi normanaturale)
‖x(ε) − x‖‖x‖ ≤ κ(A)
(‖εf‖‖b‖ +
‖εF‖‖A‖
)
Il rapporto tra l’errore relativo sui risultati e l’errore relativo sui dati risultaminore o uguale di κ(A)
SISTEMI LINEARI QUADRATI – p.11/47
Condizionamento: proprietà
κ(A) ≥ 1: infatti κ(A) = ‖A−1‖‖A‖ ≥ ‖A−1A‖ = ‖I‖ = 1
Non ci sono relazioni tra ordine del sistema e numero di condizionamento:
A1 =
1 −1
1 −1 + 10−5
,
κ∞(A1) ' 4 · 105: matrice piccola e malcondizionata
A2 =
1 −1 −1 . . . −1
0 1 −1 . . . −1
0 0 1 . . . −1
.... . .
...
0 . . . 1
, A−12 =
1 1 2 4 . . . 2n−2
0 1 1 2 . . . 2n−3
0 0 1 . . . 2n−4
.... . .
...
0 . . . 1
κ∞(A2) = ‖A2‖‖A−12 ‖ = n · 2n−1: matrice grande e malcondizionata
SISTEMI LINEARI QUADRATI – p.12/47
Condizionamento: proprietà
Non ci sono relazioni tra determinante e numero di condizionamento:det(A2) = 1≫ 0, ma κ∞(A2) = n2n−1 malcondizionata
A3 =
110
0 . . . 0
0 110
. . . 0
.... . .
...
0 . . . 110
, A−13 =
10 0 . . . 0
0 10 . . . 0
.... . .
...
0 . . . 10
,
det(A3) = 10−n, κ∞(A3) = 1: matrice quasi singolare e bencondizionata
Non è semplice calcolare κ(A): occorre calcolare A−1 che equivale a risolveren sistemi lineari ⇒ è importante stimare κ(A).
SISTEMI LINEARI QUADRATI – p.13/47
Dato Ax = b, come calcolare x?
Il metodo di Cramer è inefficiente.
Esempio: n = 30 (sistema medio-piccolo)Occorrono n + 1 determinanti di ordine n,per ciascuno occorre sommare n! termini.In ognuno ci sono n − 1 prodotti, complessivamente occorrono:(n + 1)n!(n − 1) prodotti ovvero per n = 30 circa 1034 prodotti e supponendoun tempo di 10−6 secondi per prodotto, si ottiene un tempo di circa1028sec ' 3 · 1020 anni
In qualche caso il calcolo è semplice: sistemi con matrici triangolari inferiori osuperiori.
A =
2 0 0
2 4 0
3 −2 3
:
2x = 4
2x + 4y = 8
3x − 2y + 3z = 10
x = 2
y = 8−2·24
3x − 2y + 3z = 10
x = 2
y = 1
z = 10−3·2+2·13
= 2
SISTEMI LINEARI QUADRATI – p.14/47
Sistemi triangolari
Ly = b con L ∈ Rn×n triangolare inferiore non singolare
`11y1 = b1
`21y1 + `22y2 = b2
· · ·`i1y1 + · · · + `iiyi = bi
· · ·`n1y1 + · · · · · · + `nnyn = bn
⇒
y1 = b1/`11
y2 = (b2 − `21y1)/`22
· · ·yi = [bi − (`i1y1 + · · · + `i,i−1yi−1)]/`ii
· · ·yn = (bn − ∑n−1
j=1 `njyj)/`nn
Algoritmo di Sostituzione in avanti
1. Per i = 1, 2, . . . , n
1. yi := bi
2. Per j = 1, 2, . . . , i − 1
yi := yi − `ijyj
3. yi := yi/`ii
SISTEMI LINEARI QUADRATI – p.15/47
Sistemi triangolari
Rx = y con R ∈ Rn×n triangolare superiore non singolare
r11x1 + r12x2 · · · + r1nxn = y1
r22x2 + · · · + r2nxn = y2
...
rnnxn = yn
Algoritmo di Sostituzione all’indietro
1. Per i = n, n − 1, . . . , 1
1. xi := yi
2. Per j = i + 1, i + 2, . . . , n
xi := xi − rijxj
3. xi := xi/rii
costo computazionale: si calcola il numero di prodotti/divisioni eseguite infunzione di n e si assume come costo il termine con la potenza più alta:
n∑
i=2
(i − 1) =
n−1∑
i=1
i =n(n − 1)
2prodotti
n divisioni ⇒ n2
2operazioni
SISTEMI LINEARI QUADRATI – p.16/47
Metodi diretti: fattorizzazione LRIdea Base: trasformare il problema complesso in uno equivalente più semplice.
Se esistono L ∈ Rn×n triangolare inferiore e R ∈ Rn×n triangolare superioretali che A = LR (fattorizzazione LR) allora:
AX = b ⇔ LRx = b ⇔
Ly = b
Rx = y
Esiste la fattorizzazione LR e se esiste è unica?
`11 0 · · · 0
`21 `22. . .
......
. . . 0
`n1 `n2 · · · `nn
·
r11 r12 · · · r1n
0 r22 · · · r2n
.... . .
. . ....
0 · · · 0 rnn
=
a11 · · · a1n
......
......
an1 · · · ann
n2 equazioni non lineari in n2 + n incognite ⇒ si fissano n incognite
r11 = r22 = · · · = rnn = 1: fatt. alla Crout`11 = `22 = · · · = `nn = 1: fatt. alla Dolittle
SISTEMI LINEARI QUADRATI – p.17/47
Metodi diretti: fattorizzazione LRTEOREMA di FattorizzazioneSia A ∈ Rn×n, se det(Ak) 6= 0, k = 1, 2, . . . , n − 1,(Ak minore principale di ordine k),allora esiste ed è unica la fattorizzazione alla Dolittle e det(A) =
∏ni=1 rii
(si dimostra per induzione su k)
Le ipotesi del teorema non sono di facile verifica
Sono soddisfatte dalle matrici simmetriche e definite positive e dalle matrici adiagonale dominante in senso debole e non singolari
SISTEMI LINEARI QUADRATI – p.18/47
Metodo di eliminazione di GaussCome passare da AX = b ⇔ LRx = b ? Il metodo di Gauss formalizza l’idea dirisolvere il sistema eliminando via via le incognite dalle equazioni e trasforma ilsistema assegnato in uno equivalente del tipo Rx = y.
1x + 4y + 7z = 1
2x + 5y + 8z = 1
3x + 6y + 11z = 1
→
x + 4y + 7z = 1
−3y − 6z = −1
−6y − 10z = −2
→
x + 4y + 7z = 1
−3y − 6z = −1
2z = 0
I
II
III
(1)
→
I
II = II − 21I
III = III − 31I
(2)
→
I
II
III = III − −6−3
II
(1) Elimina x dalla II[III] eq.: si moltiplica la I eq. per il coefficiente della x nellaII[III] eq. e si divide per il coefficiente di x della I eq. Si sottrae dalla II[III]
eq.
(2) Elimina y dalla III eq.: si moltiplica la II eq. per il coefficiente della y nellaIII eq. e si divide per il coefficiente di y della II eq. Si sottrae dalla III eq.
SISTEMI LINEARI QUADRATI – p.19/47
Metodo di eliminazione di GaussIn termini matriciali
1 4 7 1
2 5 8 1
3 6 11 1
→
1 4 7 1
0 −3 −6 −1
0 −6 −10 2
→
1 4 7 1
0 −3 −6 −1
0 0 2 0
[
A(0)|b(0)] (1)
→[
A(1)|b(1)] (2)
→[
A(2)|b(2)]
Anche le trasformazioni possono essere descritte in termini matriciali:
A(0) =
1 4 7
2 5 8
3 6 11
M1 :=
1 0 0
− 21
1 0
− 31
0 1
, M1A(0) = A(1) =
1 4 7
0 −3 −6
0 −6 −10
M2 :=
1 0 0
0 1 0
0 −−6−3
1
, M2A(1) = A(2) =
1 4 7
0 −3 −6
0 0 2
SISTEMI LINEARI QUADRATI – p.20/47
Metodo di eliminazione di GaussIn generale, dato AX = b si vogliono costruire M1, M2, . . . , Mn−1 tali che
Mn−1 · Mn−2 · · ·M2 · M1 · A = R
⇔ Mn−1 · Mn−2 · · ·M2 · M1 · Ax = Mn−1 · Mn−2 · · ·M2 · M1 · bRx = Mn−1 · Mn−2 · · ·M2 · M1 · b
In termini algoritmici
[A(0)|b(0)
]:= [A|b]
[A(k)|b(k)
]:= Mk
[A(k−1)|b(k−1)
], k = 1, . . . , n − 1
A(k−1) → A(k) := MkA(k−1)
0
k
k →
0k+1
k+1
SISTEMI LINEARI QUADRATI – p.21/47
Metodo di eliminazione di Gauss
A(k−1) :=
a(k−1)11 . . . . . . · · · a
(k−1)1n
0. . .
...... 0 a
(k−1)kk
......
......
. . ....
0 . . . 0 a(k−1)nk
. . . a(k−1)nn
Si vuol trovare Mk t.c.
Mk ·
a(k−1)1,k
...
a(k−1)k,k
a(k−1)k+1,k
...
a(k−1)nk
=
a(k−1)1,k
...
a(k−1)k,k
0
...
0
=:
a(k)1,k
...
a(k)k,k
0
...
0
SISTEMI LINEARI QUADRATI – p.22/47
Metodo di eliminazione di Gauss
Mk :=
1 0 . . . 0 . . . 0
0. . .
......
... 1...
... −a(k−1)k+1,k
a(k−1)k,k
1...
......
. . ....
0 . . . 0 −a(k−1)n,k
a(k−1)k,k
0 . . . 1
, α(k) :=
0
...
0
a(k−1)k+1,k
a(k−1)k,k
...a(k−1)n,k
a(k−1)k,k
⇒ Mk := I − α(k)
eTk con e
Tk = [0, . . . , 1, 0, . . . , 0]T
k
a(k−1)k,k
: elemento pivot, è cruciale per lo svolgimento del metodo:
Il metodo di Gauss si interrompe se l’elemento pivot a(k−1)k,k
è zero.
Le ipotesi del teorema di fattorizzazione (det(Ak) 6= 0, k = 1, 2, . . . , n − 1)
garantiscono che a(k−1)k,k
6= 0
SISTEMI LINEARI QUADRATI – p.23/47
Metodo di eliminazione di Gauss
A(0) =
1 4 7
2 5 8
3 6 11
, α(1) :=
0
21
31
, α(1)
eT1 =
0
2
3
[1, 0, 0] =
0 0 0
2 0 0
3 0 0
M1 = I − α(1)
eT1 =
1 0 0
−2 1 0
−3 0 1
, A(1) := M1A(0) =
1 0 0
−2 1 0
−3 0 1
1 4 7
2 5 8
3 6 11
A(1) =
1 4 7
0 −3 −6
0 −6 −10
, α(2) =
0
0
−6−3
, α(2)
eT2 =
0
0
2
[0, 1, 0] =
0 0 0
0 0 0
0 2 0
M2 = I − α(2)
eT2 =
1 0 0
0 1 0
0 −2 1
, A(2) = M2A(1) =
1 4 7
0 −3 −6
0 0 2
= R
SISTEMI LINEARI QUADRATI – p.24/47
Metodo di eliminazione di GaussIl metodo di Gauss da’ la fattorizzazione LR di A:
Mn−1 · Mn−2 · · ·M2 · M1 · A = R
A = (Mn−1 · Mn−2 · · ·M2 · M1)−1R
A = M−11 · M−1
2 · · ·M−1n−1
︸ ︷︷ ︸
·R
L := M−11 · M−1
2 · · ·M−1n−1
A = LR Fattorizzazione LR di A
M−1i
(a)
=I + α
(i)e
Ti ⇒ L
(b)
=I +
n−1∑
i=1
α(i)
eTi
eTi α
(j) = 0, se i ≤ j ⇒ (a) (I − α(i)
eTi )(I + α
(i)e
Ti ) = I − α
(i)e
Ti α
(i)e
Ti = I
(b) M−11 · M−1
2 · · ·M−1n−1 = (I + α
(1)e
T1 )(I + α
(2)e
T2 ) · · ·
= (I + α(1)
eT1 + α
(2)e
T2 + α
(1)e
T1 α
(2)e
T2 ) · · ·
= (I + α(1)
eT1 + α
(2)e
T2 )(I + α
(3)e
T3 ) · · ·
SISTEMI LINEARI QUADRATI – p.25/47
Metodo di eliminazione di Gauss
L =
1 0 . . . 0
a(0)21
a(0)11
1 0...
a(0)31
a(0)11
a(1)32
a(1)22
. . . 0 0
......
. . . 1 0
a(0)n1
a(0)11
a(1)n2
a(1)22
. . .a(n−2)n,n−1
a(n−2)n−1,n−1
1
R =
a(0)11 a
(0)12 . . . . . . a
(0)1n
0 a(1)22 . . . a
(1)2n
0 0 a(2)33 . . . a
(2)3n
.... . .
. . ....
.... . .
. . ....
0 . . . 0 a(n−1)nn
Ax = b ⇔ LRx = b ⇔
Ly = b
Rx = y⇔ Rx = L−1b
Metodo di Gauss
Mn−1 · Mn−2 · · ·M2 · M1 · A︸ ︷︷ ︸
x = Mn−1 · Mn−2 · · ·M2 · M1︸ ︷︷ ︸
·b
R x = L−1 b
Il sistema lineare Ly = b viene implicitamente risolto
SISTEMI LINEARI QUADRATI – p.26/47
Metodo di eliminazione di GaussIn pratica l’area di memoria di A viene usata per memorizzarvi R e gli elementisignificativi di L:ESEMPIO
A(0) := A =
1 −1 2
−1 5 4
2 4 14
α(1) :=
0
−1
2
,
A(1) := M1A(0) =
1 −1 2
0 4 6
0 6 10
;
1 −1 2
−1 4 6
2 6 10
α(2) =
0
0
64
,
A(2) = M2A(1) =
1 −1 2
0 4 6
0 0 1
,
1 −1 2
−1 4 6
2 64
1
SISTEMI LINEARI QUADRATI – p.27/47
Algoritmo di Gauss
1. Per k = 1, . . . , n − 1 (passi)
1. Se akk = 0 allora (controllo pivot)
→ stop
Altrimenti
→ per i = k + 1, . . . , n (righe)
1. η :=aik
akk
2. aik := η (colonna di L)
3. Per j = k + 1, . . . , n (colonne)
1. aij := aij − ηakj (riga di R)Non si costruiscono esplicitamente le matrici di trasformazione Mk
Si occupano n2 locazioni (area per memorizzare A);
Numero di operazioni (prodotti e divisioni):
n−1∑
k=1
(n − k)2 +
n−1∑
k=1
(n − k) =(n − 1)n(2n − 1)
6+
(n − 1)n
2' n3
3
SISTEMI LINEARI QUADRATI – p.28/47
Algoritmo di Gauss
0
k
ak,k(k−1)
X
X
0
k
akk(k−1)
0
0
i
aij(k−1)−η a
kj(k−1)
j
ai,j(k)
0
k
akk(k−1)
i
η=aik(k−1)/a
kk(k−1)
j
aij(k)
colonna
L
Ak−1 Ak
riga i → riga i − η riga k
↓ → riga i ↓[
0, . . . , 0, a(k−1)ik
, . . . , a(k−1)in
]
→[
0, . . . , 0, a(k−1)ik
, . . . , a(k−1)in
]
− η[
0, . . . , 0, a(k−1)kk
, . . . , a(k−1)kn
]
=[
0, . . . , 0, 0, a(k)k,k+1, . . . , a
(k)kn
]
SISTEMI LINEARI QUADRATI – p.29/47
Algoritmo di Gauss
L’algoritmo di Gauss fornisce la fattorizzazione LR di A. Per risolvere ilsistema Ax = b
occorre calcolare y = L−1b: si può aggiungere una colonna ad A ememorizzare b nella colonna n + 1 ed applicare le trasformazioni (1.1.3.1)fino alla colonna n + 1 (Per j = k + 1, . . . , n, n + 1)
rimane da risolvere Rx = y con l’algoritmo sostituzione all’indietro.
analogamente si possono risolvere in parallelo m sistemi lineari constessa matrice A
Il metodo di Gauss si arresta anche per matrici semplici:
A =
0 1
2 3
,
x2 = b1
2x1 + 3x2 = b2
applicato ad A si arresta subito, ma il sistema lineare è di immediata soluzionee il sistema equivalente ha matrice già triangolare superiore:
2x1 + 3x2 = b2
x2 = b1A =
2 3
0 1
,
SISTEMI LINEARI QUADRATI – p.30/47
Metodo di Gauss: pivoting e scaling
In generale
se a(k−1)kk
= 0 ⇒ il metodo di Gauss si interrompeIdea: cambiare di posto due righe: la k−esima con una delle righe successive
per la quale risulta a(k−1)ik
6= 0, i = k + 1, . . . , n.
(Oss: Se det(A) 6= 0 e a(k−1)kk
= 0 necessariamente qualche elemento
a(k−1)ik
, i = k + 1, . . . , n della colonna k−esima deve essere non nullo.)
Se a(k−1)kk
6= 0 ma molto piccolo (in valore assoluto) rispetto agli altri elementi⇒ l’algoritmo di Gauss può essere instabile a causa della 1.1.1
⇒ Per estendere il metodo di Gauss a matrici non singolari ed assicurare unamigliore stabilità al generico passo k−esimo si sceglie l’elemento pivot seguendouna delle due strategie:
Pivoting totale: si sceglie come elemento pivot il massimo in modulo tra tutti glielementi della matrice che rimangono da trasformare
Pivoting parziale: si sceglie come elemento pivot il massimo in modulo deglielementi di una colonna
SISTEMI LINEARI QUADRATI – p.31/47
Metodo di Gauss: pivot totale
Si sceglie a(k−1)pq con k ≤ p, q ≤ n tale che
|a(k−1)pq | = max
k≤i,j≤n|a(k−1)
ij |
se p 6= k si scambiano le righe k − p e le colonne k − q:
0
k
k
p ap,q(k−1)
q
ak,k(k−1)
∑n−1k=1 (n − k + 1)2
=n(n+1)n(2n+1)
6− 1
' n3
3confronti
Il metodo di Gauss con pivot totale risulta stabile ma in in pratica a causa dei
confronti ( n3
3) il tempo di esecuzione è circa raddoppiato
⇒ non è competitivo
SISTEMI LINEARI QUADRATI – p.32/47
Metodo di Gauss: pivot parziale
Si sceglie a(k−1)pk
con k ≤ p ≤ n tale che
|a(k−1)pk
| = maxk≤i≤n
|a(k−1)ik
|
se p 6= k si scambia la riga k−esima con la p−esima:
0
k
p
k
ak,k(k−1)
ap,k(k−1)
∑n−1k=1 (n − k + 1)
=n(n−1)
2− 1
' n2
2confronti
La stabilità del metodo di Gauss con pivot parziale non è dimostrabile masperimentalmente accertata; per n grande il tempo di esecuzione è circa lostesso di Gauss classico ⇒ risulta il metodo più usato (tra i metodi diretti persistemi lineari)
SISTEMI LINEARI QUADRATI – p.33/47
Metodo di Gauss: pivot parziale
In termini matriciali:
i < j : I = [e1, . . . , ei, . . . , ej , . . . , en], → P := [e1, . . . , ej , . . . , ei, . . . , en]
matrice identita′ matrice di permutazione
scambio di righe i − j
PA =
a11 . . . a1n
......
aj1 . . . ajn
......
ai1 . . . ain
......
an1 . . . ann
scambio di colonne i − j
AP =
a11 . . . a1j . . . a1i . . . a1n
......
......
......
......
......
......
an1 . . . anj . . . ani . . . ann
SISTEMI LINEARI QUADRATI – p.34/47
Metodo di Gauss: pivot parziale
Passo k−esimo. Si scambiano le righe k − p con Pk := [e1, . . . , ep, . . . , ek, . . . , en]:
A(k−1) → PkA(k−1)
0k
k
ak,k(k−1)
ap,k(k−1) p
0
k
ap,k(k−1)
ak,k(k−1)
Algoritmo di Gauss con pivot parziale
[A(0)|b(0)
]:= [A|b]
[A(k)|b(k)
]:= MkPk
[A(k−1)|b(k−1)
], k = 1, . . . , n − 1
SISTEMI LINEARI QUADRATI – p.35/47
Metodo di Gauss: pivot parziale
In altri terminiMn−1 · Pn−1 · · ·M1 · P1 · A = R
P := Pn−1 · · ·P1
(Mn−1 · Pn−1 · · ·M1 · P1)P−1
︸ ︷︷ ︸P · A = R
L−1 := (Mn−1 · Pn−1 · · ·M1 · P1)P−1
L := P (P1M−11 · · ·Pn−1M−1
n−1) (PiPi = I)
⇒ L−1PA = R ⇒ PA = LR Fattorizzazione LR di PA
Gauss con pivot parziale si arresta solo se:A(k−1)
k
k
0
0
⇒ det(A(k−1)) = 0
A(k−1) = Mk−1 · Pk−1 · · ·M1 · P1 · Adet(Mi) = 1, det(Pi) = ±1
det(A(k−1)) = 0 ⇔ det(A) = 0
TEOREMA
Se det(A) 6= 0 esiste una matrice di permutazione P tale che PA = LR nelsenso alla Dolittle
SISTEMI LINEARI QUADRATI – p.36/47
Metodo di Gauss: pivot parziale
Per il sistema lineare
Ax = b ⇔ PAx = Pb ⇔ LRx = Pb ⇔
Ly = Pb
Rx = y
Metodo di Gauss con Pivot parziale
Ax = b
Mn−1Pn−1 · · ·M1P1Ax = Mn−1Pn−1 · · ·M1P1b
Mn−1Pn−1 · · ·M1P1A︸ ︷︷ ︸
x = (Mn−1Pn−1 · · ·M1P1)P−1
︸ ︷︷ ︸P b
R x = L−1 P b
Il sistema lineare Ly = Pb viene implicitamente risolto
SISTEMI LINEARI QUADRATI – p.37/47
Algoritmo di Gauss con pivot parziale
1. Per k = 1, . . . , n − 1 (passi)
1. Trovare p ∈ {k, k + 1, . . . , n} tale che
|apk| = maxk≤i≤n
|aik|2. rk := p
3. Scambia akj e apj, j = k, . . . , n
4. wj := akj, j = k + 1, . . . , n
5. Per i = k + 1, . . . , n (righe)
1. η :=ai,k
akk
2. ai,k := η (colonna di L)
3. Per j = k + 1, . . . , n (colonne)
1. aij := aij − ηwj (riga di R)
L’istruzione 1.4 serve per ridurre il tempo di accesso alla memoria nella 1.5.3.1.
SISTEMI LINEARI QUADRATI – p.38/47
Algoritmo di Gauss con pivot parziale
Non si costruiscono esplicitamente le matrici di trasformazione Mk, Pk
Si occupano n2 locazioni (area per memorizzare A), le informazioni perricostruire L e P sono date da α
(i) (nelle colonne di A) e dar = [r1, . . . , rn−1]T t.c Pk = [e1, . . . , erk
, . . . , ek, . . . , en−1]
Numero di operazioni: ' n3
3
Per risolvere il sistema Ax = b occorre calcolare y = L−1Pb:si può aggiungere una colonna ad A e applicare le trasformazioni fino allacolonna n + 1 (Per j = k + 1, . . . , n, n + 1)
Occorre controllare |apk| nella 1.1.
In teoria tra 1.1 e 1.2 c’è un controllo del tipo:Se |apk| = 0 → stop
altrimenti prosegui
Non ha senso questo tipo di controllo lavorando con numeri floating-point.È preferibile:
Se |apk| ≤ tolleranza → stop
altrimenti prosegui
viene spesso usato tolleranza := εM‖A‖∞
SISTEMI LINEARI QUADRATI – p.39/47
Algoritmo di Gauss con pivot parziale
Ancora meglio in questa formaflag:= 0 (success)
...
Se |apk| = 0
flag:= −k (error)
stop
altrimenti se |apk| ≤ tolleranza
flag:= +k (warning)
prosegui
...
SISTEMI LINEARI QUADRATI – p.40/47
Metodo di fattorizzazione direttaDa LR = A, caso n = 3:
1 0 0
`421 1 0
`531 `832 1
r111 r2
12 r313
0 r622 r7
23
0 0 r933
=
a11 a12 a13
a21 a22 a23
a31 a32 a33
1) r11 = a11 4) `21r11 = a21
2) r12 = a12 1 5) `31r11 = a31 2
3) r13 = a13
6) `21r12 + r22 = a22 3 8) `31r12 + `32r22 = a32 4
7) `21r13 + r23 = a23
9) `31r13 + `32r23 + r33 = a32 5
I termini noti si usano una sola voltain corrispondenza degli elementiche si costruiscono
(5)
(2)
(4)
(3)
(1)
SISTEMI LINEARI QUADRATI – p.41/47
Metodo di fattorizzazione direttaDa LR = A, caso n = 3:
1 0 0
`421 1 0
`531 `832 1
r111 r2
12 r313
0 r622 r7
23
0 0 r933
=
a11 a12 a13
a21 a22 a23
a31 a32 a33
1) r11 = a11 4) `21r11 = a21
2) r12 = a12 1 5) `31r11 = a31 2
3) r13 = a13
6) `21r12 + r22 = a22 3 8) `31r12 + `32r22 = a32 4
7) `21r13 + r23 = a23
9) `31r13 + `32r23 + r33 = a32 5
I termini noti si usano una sola voltain corrispondenza degli elementiche si costruiscono
(5)
(2)
(4)
(3)
(1)
SISTEMI LINEARI QUADRATI – p.41/47
Metodo di fattorizzazione direttaIn generale → Pavimentazione alla Dolittle
(2) (4)
01
1
11 0
(3)
(1)
(2) (4)
(3)
(1)
L R A
Richiede n3
3operazioni e circa n2 locazioni
In pratica è instabile e richiede strategie di pivoting
SISTEMI LINEARI QUADRATI – p.42/47
Metodo di fattorizzazione direttaIn generale → Pavimentazione alla Dolittle
(2) (4)
01
1
11 0
(3)
(1)
(2) (4)
(3)
(1)
L R A
Richiede n3
3operazioni e circa n2 locazioni
In pratica è instabile e richiede strategie di pivoting
SISTEMI LINEARI QUADRATI – p.42/47
Metodo di fattorizzazione diretta
Per A simmetrica AT = A: siano
D :=
r11 0 . . . 0
0 r22
. . ....
.... . .
. . . 0
0 . . . 0 rnn
U :=
1 r12r11
r13r11
. . . r1n
r11
0 1 r23r22
. . . r2n
r22
.... . .
...
0 . . . 0 1
R = DU, A = LDU e AT = UT DLT ;
ma A = AT ⇒ L = UT , `ij =
uji =rji
rjj, i > j
1, i = j
0, i < jIl calcolo di L risulta semplificato:
1 0 0
`2′
21 1 0
`3′
31 `5′
32 1
r111 r2
12 r313
0 r422 r5
23
0 0 r633
=
a11 a12 a13
a21 a22 a23
a31 a32 a33
Richiede circa n3
6operazioni
SISTEMI LINEARI QUADRATI – p.43/47
Metodo di fattorizzazione diretta
Per A simmetrica AT = A: siano
D :=
r11 0 . . . 0
0 r22
. . ....
.... . .
. . . 0
0 . . . 0 rnn
U :=
1 r12r11
r13r11
. . . r1n
r11
0 1 r23r22
. . . r2n
r22
.... . .
...
0 . . . 0 1
R = DU, A = LDU e AT = UT DLT ;
ma A = AT ⇒ L = UT , `ij =
uji =rji
rjj, i > j
1, i = j
0, i < jIl calcolo di L risulta semplificato:
1 0 0
`2′
21 1 0
`3′
31 `5′
32 1
r111 r2
12 r313
0 r422 r5
23
0 0 r633
=
a11 a12 a13
a21 a22 a23
a31 a32 a33
Richiede circa n3
6operazioni
SISTEMI LINEARI QUADRATI – p.43/47
Fattorizzazione di Choleski
Per A simmetrica e definita positiva AT = A e xT Ax > 0 ∀x 6= 0 x ∈ Rn:TEOREMA
Sia A ∈ Rn×n simmetrica e definita positiva. Allora det (Ak) > 0, k = 1, . . . , n
Dal Terorema di fattorizzazione det (A) =∏n
i=1 rii ⇒ che rii > 0, si pone
D12 :=
√r11 0 . . . 0
0√
r22
. . ....
.... . .
. . . 0
0 . . . 0√
rnn
= diag(√
r11, . . . ,√
rnn)
A = LDU = LDLT = LD12 D
12 LT , posto S := LD
12 ⇒ A = SST
Da notare che se M = NNT , det (N) 6= 0 allora M è simmetrica e definita positiva:∀x ∈ Rn, x 6= 0 : xT Mx = xT NNT x = ‖NT x‖2
2 > 0 ⇒TEOREMA di Choleski
A ∈ Rn×n è simmetrica e definita positiva se e solo se ammette unafattorizzazione del tipo A = SST con S triangolare inferiore.Fattorizzazione di Choleski
SISTEMI LINEARI QUADRATI – p.44/47
Fattorizzazione di Choleski
Per A simmetrica e definita positiva AT = A e xT Ax > 0 ∀x 6= 0 x ∈ Rn:TEOREMA
Sia A ∈ Rn×n simmetrica e definita positiva. Allora det (Ak) > 0, k = 1, . . . , n
Dal Terorema di fattorizzazione det (A) =∏n
i=1 rii ⇒ che rii > 0, si pone
D12 :=
√r11 0 . . . 0
0√
r22
. . ....
.... . .
. . . 0
0 . . . 0√
rnn
= diag(√
r11, . . . ,√
rnn)
A = LDU = LDLT = LD12 D
12 LT , posto S := LD
12 ⇒ A = SST
Da notare che se M = NNT , det (N) 6= 0 allora M è simmetrica e definita positiva:∀x ∈ Rn, x 6= 0 : xT Mx = xT NNT x = ‖NT x‖2
2 > 0 ⇒TEOREMA di Choleski
A ∈ Rn×n è simmetrica e definita positiva se e solo se ammette unafattorizzazione del tipo A = SST con S triangolare inferiore.Fattorizzazione di Choleski
SISTEMI LINEARI QUADRATI – p.45/47
Fattorizzazione di Choleski
s11 0 · · · 0
s21 s22
. . ....
.... . . 0
sn1 sn2 · · · snn
·
s11 s21 · · · sn1
0 s22 · · · sn2
.... . .
. . ....
0 · · · 0 snn
=
a11 · · · a1n
......
......
an1 · · · ann
aii =
i∑
k=1
s2ik ⇒ sii =
√√√√aii −
i−1∑
k=1
s2ik
aij =
j∑
k=1
siksjk, i > j ⇒ sij =
aij −j−1∑
k=1
siksjk
/sjj
Algoritmo di Choleski (per colonne)1. Per j = 1, 2, . . . , n
1. sjj :=√
ajj − ∑j−1k=1 s2
jk
2. Per i = j + 1, . . . , n
sij = (aij −j−1∑
k=1
siksjk)/sjj
SISTEMI LINEARI QUADRATI – p.46/47
Fattorizzazione di Choleski
s11 0 · · · 0
s21 s22
. . ....
.... . . 0
sn1 sn2 · · · snn
·
s11 s21 · · · sn1
0 s22 · · · sn2
.... . .
. . ....
0 · · · 0 snn
=
a11 · · · a1n
......
......
an1 · · · ann
aii =
i∑
k=1
s2ik ⇒ sii =
√√√√aii −
i−1∑
k=1
s2ik
aij =
j∑
k=1
siksjk, i > j ⇒ sij =
aij −j−1∑
k=1
siksjk
/sjj
Algoritmo di Choleski (per colonne)1. Per j = 1, 2, . . . , n
1. ajj :=√
ajj − ∑j−1k=1 ajk
2
2. Per i = j + 1, . . . , n S viene memorizzata nella parte inferiore di A
aij = (aij −j−1∑
k=1
aikajk)/ajj
SISTEMI LINEARI QUADRATI – p.46/47
Fattorizzazione di Choleski
n2 locazioni per S senza perdere A (basta memorizzare la diagonale di A inun vettore ausiliario)
Richiede n3
6operazioni
L’algoritmo di Choleski è stabile
Occorre controllare se radicando > 0 nella 1.1. Si fa un controllo del tiporadicando ≥ tolleranza
per evitare che nella 1.2.1 si divida per ajj ∼ 0
L’algoritmo di Choleski è usato per verificare se A simmetrica è definita positiva
SISTEMI LINEARI QUADRATI – p.47/47
Fattorizzazione di Choleski
n2 locazioni per S senza perdere A (basta memorizzare la diagonale di A inun vettore ausiliario)
Richiede n3
6operazioni
L’algoritmo di Choleski è stabile
Occorre controllare se radicando > 0 nella 1.1. Si fa un controllo del tiporadicando ≥ tolleranza
per evitare che nella 1.2.1 si divida per ajj ∼ 0
L’algoritmo di Choleski è usato per verificare se A simmetrica è definita positiva
SISTEMI LINEARI QUADRATI – p.47/47