1. No category

Elementære eliminasjonsmatriser
Gitt en vektor a = [a1 , . . . , an ]T ,

1 ···
0
0
 .. . .
..
..
 .
.
.
.

 0 ···
1
0
Mk = 
 0 · · · − ak+1
1
ak

..
..
 .. . .
.
 .
.
.
0 ···
− aank
en matrise

··· 0
. . .. 
. . 

··· 0 
,
··· 0 

. . .. 
. . 
0 ··· 1
1/17
 
a1
 ..  
 .  
 

 a  
Mk  k  = 
 ak+1  
 .  
 ..  
an

ak 6= 0,
kalles elementære eliminasjonsmatriser
eller Gauss transformasjon.
ak kalles
Fakta om Mk :
• Mk er nedre triangulær med 1-ere p˚
a diagonalen, derfor ikke-singulær
• Mk = I − meTk , m = [0, . . . , 0, mk+1 , . . . mn ]T
• Mk−1 = I + meTk = Lk

a1
.. 
. 

ak 

0 
.. 
. 
0
pivot.
JJ
II
J
I
Back
Close
• Hvis j > k, Mj = I − ueTj , da:
Mk Mj = I − meTk − ueTj + meTk ueTj = I − meTk − ueTj ,
P
siden eTk u = eTk nl=j+1 ul el = 0.
2/17
Det samme for Lk Lj , siden
Lj = Mj−1 = I + meTj .
Merk rekkefølge!
JJ
II
J
I
Back
Close
Gauss eliminasjon, LU faktorisering
3/17
Ax = b
• Multiplisere begge sider med M1 , med a1,1 som pivot, slik at [a1,1 , a2,1 , . . . , an,1 ]T 7→ [a1,1 , 0, . . . , 0]T .
• Mult. med M2 s˚
a at elementene under diagonalen er sett til null i kolonne 2.
• ...
Etter n − 1 skritt, vi har
Mn−1 Mn−2 · · · M1 Ax = Mn−1 Mn−2 · · · M1 b
hvor
U = M A = Mn−1 Mn−2 · · · M1 A
er øvre triangulær.
Vi kan finne x ved baklengs substitusjon i U x = M b.
Denne prosessen kalles
Gauss eliminasjon
JJ
II
J
I
Back
Close
Eksempel
Regn ut Gauss eliminasjons metode for
2x1 + 3x2 + 4x3 + 5x4
x1 + x2 − x3
3x1 − x2 + 3x3 − x4
x2 − 3x3 + x4
4/17
=
=
=
=
−5
0
0
0.
JJ
II
J
I
Back
Close
Hvorfor kalles dette ogs˚
a LU ?
5/17
Vi har sett at
M A = U ⇒ A = M −1 U = LU
M = Mn−1 · · · M1 er nedre triangulær, s˚
a inversen er ogs˚
a n.t.
−1
L = M −1 = (Mn−1 · · · M1 )−1 = M1−1 · · · Mn−1
= L1 · · · Ln−1 .
Hvis vi antar at A = LU er gitt,
• først setter vi y = U x og deretter løser vi Ly = b (forlengs substitusjon)
• Vi løser U x = y (baklengs substitusjon)
Merk at den midlertidig vektor y = L−1 b = M b.
Gauss eliminasjon og LU er to sider av den samme medalje.
• LU faktorisering kan implementeres uten ˚
a modifisere b
• LU faktorisering er anbefalt n˚
ar man skal løse mange likningssystemer med samme A og
forskjellige b.
JJ
II
J
I
Back
Close
Litt om implementasjon
• De superdiagonale elementer av U erstatter A sine elementer
6/17
• De underdiagonale elementer av A (som blir null) brukes til ˚
a lagre L sine elementer
Denne prosedyren kalles
factorization in place (faktorisering p˚a plass).
Algoritme: LU faktorisering m/ Gauss
eliminasjon
Algoritme: p˚
a plass LU faktorisering m/
Gauss eliminasjon
for k = 1 to n − 1
if ak,k = 0, stop
for i = k + 1 to n
li,k = ai,k /ak,k
end
for j = k + 1 to n
for i = k + 1 to n
ai,j = ai,j − li,k ak,j
end
end
end
for k = 1 to n − 1
if ak,k = 0, stop
for i = k + 1 to n
ai,k = ai,k /ak,k
end
for j = k + 1 to n
for i = k + 1 to n
ai,j = ai,j − ai,k ak,j
end
end
end
JJ
II
J
I
Back
Close
Pivotering
7/17
Hvis pivoten er null, da kan ikke Gauss eliminasjon utføres
ak,k = 0 ⇒ mk+1 =
an,k
ak+1,k
, . . . , mn =
er ikke definert
ak,k
ak,k
og faktorisering m˚
a stoppes (selv om A er ikke-singulær)
Et annet tilfelle er hvis pivoten ak,k er veldig liten, ak,k ≈
mi =
ai,k
,
ak,k
i = k + 1, . . . , n,
1
kan være veldig store
og hvis de andre ai,j er av moderat størrelse, vi kan tape mange signifikante siffer
Eksempler
A=
0 1
1 0
,
A=
1
1 1
,
JJ
II
J
I
Back
Close
Partial Pivoting
I prinsippet: store pivoter ⇒ sm˚
a mk og derfor mindre feil
Vi søker etter den største verdien under diagonalen i kolonne k (dermed ‘kolonnevis pivoting’).
Hvis denne er i rekke p, bytter vi ut rekker k og p og faktoriserer som vanlig.
Merk at n˚
a er
|mk | ≤ 1
8/17
Husk: bytting av rekker ⇔ permutasjoner
M A = U,
M = Mn−1 Pn−1 · · · M1 P1
hver elementære eliminasjonsmatrise etterfølger en permutasjons matrise.
La oss skrive
P = Pn−1 · · · P1
Gauss eliminasjonen m/ partial pivoting er ekvivalent til den LU faktorisering av P A:
P A = LU
L, U nedre/øvre triangulære.
Ax = b
Obs. P er ikke kjent p˚
a forhand.
Ly = P b,
U x = y.
JJ
II
J
I
Back
Close
Algoritme: p˚
a plass LU faktorisering med Gauss eliminasjon og kolonnevis pivoting
for k = 1 to n − 1
find index p s.t. |ap,k | ≥ |ai,k |, for k ≤ i ≤ n
if p 6= k, bytt ut rekker p og k
if ak,k = 0 continue with next k
% hopper over denne kolonne, alle elementer er 0 allrede
for i = k + 1 to n
ai,k = ai,k /ak,k
end
for j = k + 1 to n
for i = k + 1 to n
ai,j = ai,j − ai,k ak,j
end
end
end
9/17
JJ
II
J
I
Back
Close
Eksempel
Regn ut LU m/ kolonne pivotering for
2x1 + 3x2 + 4x3 + 5x4
x1 + x2 − x3
3x1 − x2 + 3x3 − x4
x2 − 3x3 + x4
10/17
=
=
=
=
−5
0
0
0.
JJ
II
J
I
Back
Close
Total pivoting
Kolonne pivotering er en partiell pivotering, siden pivoten er søket i en kolonne per gang.
Det finnes en annen type pivoting som
11/17
• søker det største elementet i hele submatrisen som skal prosesseres
• bytter ut kolonner og rekker slik at det største elementet kommer i pivotal plass
Denne prosedyren kalles for
total pivotering og er ekvivalent med ˚a faktorisere
P AQ = LU,
hvor P, Q er permutasjons matriser og L, U er nedre og øvre triangulære matriser.
Løsning:
Ly = P b
Uz = y
x = Qz
JJ
II
J
I
Back
Close
Stabilitet av total pivotering er i prinsippet bedre enn kolonne pivotering, men total pivotering er
noe dyrere.
I praksis, i de fleste tilfeller er kolonne pivotering god nok, dermed er den pivotering strategien som
brukes mest.
12/17
Det finnes ogs˚
a noe matriser som trenger ikke pivotering for at Gauss eliminasjons metode er stabil.
Disse er:
P
• diagonal dominante matriser: |aj,j | > i6=j |ai,j |, for j = 1, . . . , n,
• symmetrisk positiv definite matriser: A = AT og uT Au > 0 for all u 6= 0.
Selv om pivotering kan brukes p˚
a disse matriser, blir det neppe noe rekke-bytte. Med ˚
a ikke-pivotere
vi sparer en del beregningstid.
JJ
II
J
I
Back
Close
Om kompleksitet
Hvor mye koster det ˚
a løse Ax = b?
13/17
3
3
• LU faktorisering koster: ≈ n /3 flyttall multiplikasjoner, ≈ n /3 addisjoner , for en total av
2 3
n
3
• hver triangulært system koster cirka n2 multiplikasjoner og n2 addisjoner, for en total av 2n2
Det er klart at for store n, den LU-faktorisering kostnade er størst.
Merk at en matrise inversjon koster ≈ 2n3 , tre ganger mer enn LU faktorisering.
Derfor lønner det seg ˚
a regne ut produkter som
A−1 B
ved
• LU faktorisere A
• løse n forlengs/baklengs systemer (ta som b en kolonne av B om gangen).
JJ
II
J
I
Back
Close
Spesielle systemer
Opp til n˚
a har vi regnet med at A er en vilk˚
alig full (‘dense’) matrise.
Men det finnes mange matriser som har spesielle egenskaper: disse egenskaper kan brukes til ˚
a
minske beregningstid og lagring av data.
Noen spesielle tilfeller er:
14/17
• A symmetrisk, A = AT (eller ai,j = aj,i )
• A positiv definitt, xT Ax > 0, x 6= 0
• A har b˚
and β: |ai,j | = 0 for |i − j| > β,
Tridiagonale matriser har b˚
and β = 1
• A glissen (’sparse’), de fleste elementer av A er lik 0.
Symmetriske matriser kan faktoriseres med en variant av Gauss eliminasjonsmetode (Choleski faktorisering, A = LLT )
Likningssystemer med sparse matriser løses med
• direkte algoritmer som tar hensyn til 0-datastrukturen (man vil helst unng˚
a fill-in)
• iterative metoder
JJ
II
J
I
Back
Close
Symmetriske positive definite
Hvis A er symmetrisk og positiv definitt, kan man modifisere LU algoritme slik at U = LT og
dermed A = LLT .
15/17
Algoritme: p˚
a plass Choleski faktorisering
for k = 1 to n
√
ak,k = ak,k
for i = k + 1 to n
ai,k = ai,k /ak,k
end
for j = k + 1 to n
for i = k + 1 to n
ai,j = ai,j − ai,k aj,k
end
end
end
JJ
II
J
I
Back
Close
Denne kalles
Choleski faktorisering og har en del fordeler i forhold til vanlig LU :
• kvadrat-rotene er vel definerte
• ingen pivotisering
• mindre lagringsplass (lagre bare L)
• koster halvparten av LU,
1 3
n
3
16/17
(add+mult).
En variasjon: LDLT
2
• diagonale elementene av D settes til li,i
• diagonale elementene av L settes til 1
JJ
II
J
I
Back
Close
B˚
and systemer
LU faktorisering for b˚
and matriser er ikke veldig mye anneledes enn for dense matriser.
17/17
• pass p˚
a nedre og øvre grensene for loop indeksene.
Hvis man har pivotering p˚
a grunn av numerisk stabilitet, den opprinnelig b˚
and β kan ikke bli større
2β.
Generelt, b˚
and systemer trenger bare O(βn) lagringsplass og O(βn2 ) beregningsarbeide, som kan
være betidelig mindre enn O(n3 ) for β n
JJ
II
J
I
Back
Close