Examen Numerieke Analyse I, 1 lic Info, 5 juni 2000
deel 1
1. Laten u = (u1,· · · , un)T en v = (v1,· · · , vn)T twee re¨ele vektoren zijn en S := uTv het inprodukt ervan, berekend met de algoritme
S := 0; for k := 1 to n do S := S + uk∗ vk.
Laat d.m.v. een afrondfoutenanalyse zien dat voor de berekende waarde van S (in floating point met machineprecisie η) geldt: | S − fl(S) | ≤ (n + 1) η k u k2k v k2.
antwoord:
Afrondfout in de berekende waarde van het inprodukt S :=
Xn i=1
xiyi berekend met algoritme: S := 0;
for i := 1 to n do S := S + xi ∗ yi
Voor de berekende waarde van S vinden we getallen ξien ηi met | ξi|, | εi| ≤ η, i = 1 · · · n : fl(S) = x1y1(1 + ξ1) (1 + ε2) · · · (1 + εn)
+ x2y2(1 + ξ2) (1 + ε2) · · · (1 + εn) + · · ·
+ xn−2yn−2(1 + ξn−2) (1 + εn−2) · · · (1 + εn) + xn−1yn−1(1 + ξn−1) (1 + εn−1) (1 + εn) + xnyn(1 + ξn) (1 + εn)
zodat
S − fl(S) = Xn
i=1
xiyiζi
met
ζi := 1 − (1 + ξi) (1 + εi) · · · (1 + εn) en | ζi| ≤ (n − i + 2) η als nη≤ 0.1 . Bijgevolg geldt voor de voorwaartse fout:
|S − fl(S)
S | ≤ (n + 1) η
| S |
Xn i=1
| xiyi| ≤ (n + 1) η k x k2k y k2
| xTy| 2. Laat A∈ IRn×n een inverteerbare matrix zijn.
a. Als je a priori weet, dat een ontbinding van A in een product van een linksondermatrix L met diagonaalelementen L(i, i) = 1 en een rechtsbovenmatrix U mogelijk is, laat dan zien dat het mogelijk is de niettriviale elementen van L en U te berekenen in de volgorde van de Doolittle-algoritme:
for k = 1 : n do bereken U (1 : k, k) ; bereken L(k + 1 : n, k) ; end
(bereken in de k-de slag de k-de kolommen van L en U ) en schrijf die algoritme uit.
Antwoord: Tussen de elementen van A, L en U geldt, ongeacht de manier waarop ze berekend zijn, de relatie
A(i, k) =
min(i,k)X
j=1
L(i, j)U (j, k) = ( als i≤ k ) U(i, k) +iX−1
j=1
L(i, j)U (j, k) . (1)
Uit de tweede formule van (1) zien we, dat we U (i, k) voor i≤ k kunnen berekenen met
U (i, k) = A(i, k)−kX−1
j=1
L(i, j)U (j, k), (2)
als U (1 : i− 1, k) en L(1 : i − 1, i − 1) bekend zijn. Uit de eerste formule van (1) zien we dat we L(i, k) voor i > k kunnen berekenen met
L(i, k) = A(i, k)−Pij=1−1 L(i, j)U (j, k)
U (k, k) , (3)
als U (1 : k, k) en L(1 : k− 1, j) bekend zijn. De algoritme luidt dus:
for k=1:n, for i=1:k,
U(i,k) = A(i,k)-L(i,1:i-1)*U(1:i-1,k);
end
L(k+1:n,k) = A(k+1:n,k)-L(k+1:n,1:k-1)*U(1:k-1,k)
U(k, k) ;
end
(4)
b. Laat zien, dat je ook bij deze volgorde de rijverwisselingsstrategie van de Gausseliminatie kunt inbouwen (voor matrices waarbij rijverwisselingen wel nodig zijn voor de stabiliteit) en geef de aangepaste algoritme.
Antwoord: Als we de hierboven afgeleide algoritme voor de berekening van L en U vergelijken met de berekening van L en U via Gausseliminatie (zonder rijverwisseling),
L=eye(size(A));
for k=1:n-1,
L(k+1:n,k) = A(k+1:n,k)/A(k,k);
A(k+1:n,k:n) = A(k+1:n,k:n)-L(k+1:n,k)*A(k,k:n);
end
(5)
dan zien we dat de teller van (3) in het begin van de k-de slag precies de elementen van de k-de kolom van de (in k − 1 Gauss-slagen getransformeerde) matrix A bevat. Om de rijverwisselingsstrategie van de Gausseliminatie in te passen, berekenen we dus eerst deze tellers, zoeken er het element met maximale absolute waarden in en verwisselen de
bijbehorende rijen. Dit geeft dan de algoritme:
p=1:n; {array voor permutatie-indices}
for k=1:n,
for i=1:k-1,
U(i,k) = A(i,k)-L(i,1:i-1)*U(1:i-1,k);
end
Z = A(k:n,k)-L(k:n,1:k-1)*U(1:k-1,k)
{overeenkomstige kolom in de k-slag v/d Gausseliminatie}
[m,pp] = max(abs(Z)); p(k) = k-1+pp;
U(k,k) = Z(pp); Z(pp) = Z(1);
L(k+1:n) = Z(2:n-k+1)/U(k,k);
H = L(k,1:k-1); L(k,1:k-1) = L(p(k),1:k-1); L(p(k),1:k-1) = H;
H = A(k,k+1:n); A(k,k+1:n) = A(p(k),k+1:n); A(p(k),k+1:n) = H;
end
(6)
Hierin is L(i,1:i-1)*U(1:i-1,k) een inproduct van twee vectoren met lengte i-1 en A(k:n,1:k-1)*A(1:k-1,k) zijn n-k+1 inproducten .
Als je de matrix overschrijft met de factoren L en U , wordt de algoritme compacter:
p=1:n;
for k=1:n,
for i=1:k-1,
A(i,k) = A(i,k)-A(i,1:i-1)*A(1:i-1,k);
end
A(k:n,k) = A(k:n,k)-A(k:n,1:k-1)*A(1:k-1,k) [m,pp] = max(abs(A(k:n,k))); p(k) = k-1+pp;
if p(k)>k , H = A(k,:); A(k,:) = A(p(k),:); A(p(k),:) = H; end A(k,k+1:n) = A(k,k+1:n)/A(k,k);
end
(7)
Bij implementatie in Matlab moet je er wel aan denken, dat een statement van de vorm U(i,k) = A(i,k)-L(i,1:i-1)*U(1:i-1,k) crasht als i = 1 , omdat L(i,1:i-1) dan een lege matrix is.
c. Wat is het voordeel van deze algoritme t.o.v. Gausseliminatie bij gebruik op je PC, waar je re¨ele getallen in het geheugen 64 bits beslaan en waarin je floating point unit in de CPU met 96-bits registers rekent? (gebruik het resultaat van vraag 1).
Antwoord Bij een berekening van de vorm
U (i, k) = A(i, k)− L(i, 1 : i − 1) ∗ U(1 : i − 1, k);
wordt het betreffende element van U in een maal uitgerekend door middel van een in- product. Bij een goede implementatie kunnen de tussenresultaten dus in de registers blijven, zodat het inproduct dus met 96-bits precisie berekend wordt en pas daarna wordt afgerond naar 64 bits. De absolute afrondfout in het resultaat is dus begrensd door η1kLk kUk + η2(n + 1)kLk kUk met η2 ¿ η1.
3. Bewijs de interpolatieformule van Lagrange voor het geval n = 2:
Als f ∈ C2[a, b] en als x1 < x2 twee verschillende punten in het interval [a, b] zijn, dan is er bij iedere x∈ [a, b] een ξ ∈ [a, b] zodat
f (x) = f (x1) x− x2
x1− x2
+ f (x2)x− x1
x2− x1
+ 1
2(x− x1)(x− x2)f00(ξ) . (8)
Bewijs vervolgens de integratieformule en restterm, d.w.z. bewijs dat er een η ∈ [a, b] bestaat
zodat: Z
b
a f (x)dx = 12(b− a)(f(a) + f(b)) − 121(b− a)3f00(η) (9) Antwoord: Het interpolatiepolynoom van graad 1 op de interpolatiepunten x1 en x2 is
π(x) := f (x1)x− x2
x1− x2
+ f (x2)x− x1
x2− x1
.
We defini¨eren de functies g en Ry door
f (x)− π(x) − (x − x1)(x− x2)g(x) = 0 en Ry(x) := f (x)− π(x) − (x − x1)(x− x2)g(y) . De functie Ry heeft (minstens) 3 nulpunten, namelijk x1, x2 en y. De afgeleide ervan heeft (volgens Rolle) dus minstens 2 en de tweede afgeleide minstens ´e´en nulpunt tussen x1, x2
en y; noem dit nulpunt van de tweede afgeleide ξy (het hangt immers van y af). Dan geldt 0 = R00y(ξy) = f00(ξy)− 2g(y) .
Hieruit volgt (8).
Om (9) te bewijzen integreren we (8) met x1 = a en x2 = b . De integraal van π geeft
1
2(b− a)(f(a) + f(b)) (oppervlak van het trapezium onder de rechte π tussen de verticale lijnen x = a en x = b . Omdat (x − a)(x − b) (semi-)definiet is op [a, b] mogen we de middelwaardestelling van de integraalrekening toepassen en is er een η ∈ (a, b) zodat
Z b
a(x− a)(x − b)f00(ξx) dx = f00(η)
Z b
a (x− a)(x − b) dx = − 1
12(b− a)3f00(η)
4. Een Givensrotatie J(k, l, ϑ) ∈ IRm×m met l > k is een lineaire transformatie, die de vektoren in het vlak opgespannen door ek en el roteert over een hoek ϑ:
J(k, l, ϑ)(αek + βel) = (αc + βs) ek + (−αs + βc) el met c := cos ϑ en s := sin ϑ , J(k, l, ϑ) x = x , ∀x ∈ vect{ek, el}⊥.
a. schrijf de matrixrepresentatie op van J(k, l, ϑ) voor m = 5, k = 2, l = 5.
Antwoord:
J(k, l, ϑ) =
1 0 0 0 0
0 c 0 0 s
0 0 1 0 0
0 0 0 1 0
0 −s 0 0 c
met c = cos(ϑ) en s = sin(ϑ) .
b. Laten c en s zo zijn, dat voor de `-de component van een zekere vector a ∈ IRm geldt:
(J(k, `, ϑ) a)` = 0 .
Laat zien, dat deze s en c (met c ≥ 0 en dus met −12π ≤ ϑ ≤ 12π) alsvolgt berekend kunnen worden:
if| ak| > | al| then t := al/ak; c := 1/√
1 + t2; s := t ∗ c else
t := ak/al; s := sign(t)/√
1 + t2; c := | t ∗ s | end
(10)
Antwoord: (J(k, `, ϑ) a)` = − s ak+ c a` = 0 . Kwadrateren geeft c2a2` = s2a2k = (1− c2) a2k zodat c2 = a2k
a2k+ a2` en s = c a`
ak
. (11)
Aangezien we de cosinus (= c) positief willen hebben nemen we in (11) de positieve wortel. Uitdelen van ak of a` geeft dan
c = 1
q
1 + a2`/a2k, s = c a`
ak
resp. c = q| ak/a`|
1 + a2k/a2` , s = a`
ak
| ak/a`|
q
1 + a2k/a2` (12) Door eerst t te berekenen, zoals in (10), en vervolges c en s (resp. s en c) minimaliseren we het aantal benodigde flops. Bovendien minimaliseren we de op te bouwen afrondfout door altijd de wortel√
1 + t2 te berekenen met −1 ≤ t ≤ 1.
Het is duidelijk dat de rotatiehoek ϑ nooit berekend hoeft te worden.
c. We kunnen het element alk(l 6= k) van een matrix A ∈ IRm×n op nul transformeren, door de matrix van links te vermenigvuldigen met een geschikte Givensrotatie J(k, l, ϑ).
Als we daarna de resulterende matrix van links vermenigvuldigen met een Givensrotatie om een ander element nul te maken, kan het zijn dat het element alk weer een waarde ongelijk aan nul krijgt (geef een voorbeeld).
Geef aan, hoe de m− 1 elementen a2,1, a3,1, ..., am,1 van de eerste kolom van de matrix A d.m.v. een serie van m − 1 Givensrotaties nul gemaakt kunnen worden. Laat ook zien, dat een QR-ontbinding van A ∈ IRm×n gemaakt kan worden met n(n− 1)/2 + n(m − n) Givensrotaties (geef de algoritme) en bereken het aantal benodigde flops.
Antwoord:
Voorbeeld; als we het tweede en derde element van de (kolom-)vector a := (1 , 1 , 1)T ∈ IR3 op nul willen roteren, zouden we eerst de rotatie J(1, 2, ϑ) kunnen nemen met een ϑ zo dat de tweede component van b := J(1, 2, ϑ) a nul is (cos(ϑ) = 12√
2 in dat geval).
Als we vervolgens de derde component van b nul willen maken met een rotatie J(2, 3, ϕ), dan kan de tweede component van het resultaat niet nul blijven (omdat de som van de kwadraten van beide elementen niet verandert, de transformatie is orthogonaal). Als we deze nul op de tweede plaats willen bewaren, moeten we een rotatie van de vorm J(1, 3, ϕ) kiezen.
In het algemeen, als we alle componenten van een vector a ∈ IRn behalve de eerste op nul willen roteren, kunnen we vele volgordes kiezen. Veel gebruikt zijn:
for k = 2 : n, maak ak nul door rotatie J (1, k, ϕ) met geschikte hoek ϕ end for k = n :−1 : 2, maak ak nul door rotatie J (k− 1, k, ϕ) met geschikte hoek ϕ end
Voor een matrix A∈ IRm×n met m ≥ n kunnen we de vogende orde gebruiken:
for k = 1 : n,
for ` = m :−1 : k + 1,
maak A`,k nul door rotatie J (k, `, ϕ) met geschikt gekozen ϕ end
end Telling levert
Xn k=1
m− k = 1
2n(m− 1 + m − n) = nm − 1
2n(n + 1) = n(n− 1)/2 + n(m − n) Givens rotaties voor deze QR-ontbinding.
Het aantal flops per rotatie in de k-de slag is 6(n− k + 1) en dus totaal
Xn k=1
6(m− k)(n − k + 1) = n (n + 1) (3m − n − 2) .
d. Laat zien dat je met de gevonden QR-ontbinding het kleinste-kwadratenprobleem arg min
x∈IRnkAx − bk2
kunt oplossen, als de rang van A gelijk is aan n en m > n (geef ook aan hoe).
Antwoord: De bovenstaande algoritme geeft een ontbinding van A = Q R in een product van een orthogonale m× m matrix Q (met QT Q = I = Q QT) en een bovendriehoeks- matrix R ∈ IRm×n, die gelijke rang heeft als A. Deze bovendriehoeksmatrix kunnen we partitioneren in een vierkante bovendriehoeksmatrix R1 bestaande uit de eerste n rijen van R en een (m− n) × n bestaande uit nullen. Aangezien de Euclidische norm invariant is onder orthogonale transformaties geldt
kAx − bk22 =kQ R x − bk2 =kR x − QT bk2 =kµR1x 0
¶
−µd r
¶
k2 ≥ krk2 (13)
met µ
d r
¶
= QT b , d∈ IRn en r∈ IRm−n
een partitie van QT b . De eerste n componenten van de vector tussen de meest rechtse normstrepen in (13) kunnen nul gemaakt worden door de (unieke) keuze x = R−11 d of beter gezegd door het oplossen van x uit R1x = d . Hetgeen overblijft is krk2. Geen enkele keuze van x kan dit deel van het rechterlid wegwerken; voor iedere andere keuze van x dan bovenstaande blijft er een residu, dat groter is. De gevonden x geeft dus het minimum en dus de oplossing van het kleinste-kwadratenprobleem.
Als je de “economy size” QR-ontbinding A = Q Re 1 met Qe ∈ IRm×n en R1 ∈ IRn×n (als boven) neemt, zoals je die bijvoorbeeld uit de (Modified) Gram-Schmidt algoritme krijgt, en je definieert y := R1x , dan is het kleinste-kwadratenprobleem
arg min
x∈IRnkAx − bk2 = arg min
y∈IRnkQye − bk2
equivalent met de normaalvergelijkingenQeT Q y =e QeT b en dus y = QeTb (= d) .
