Numerieke Lineaire Algebra

(1)

P. de Groen

In deze cursusnota’s bij de cursus Numerieke Lineaire Algebra van het tweede jaar bachelor in de wiskunde worden de standaard algoritmen behandeld voor het oplossen van een stelsel lineaire vergelijkingen en een lineair kleinste-kwadratenprobleem. Voor een goed begrip van de werking van de algoritmen in de praktijk wordt dit voorafgegaan door een inleiding in de analyse van afrondfouten en een inleiding in Matlab als eenvoudig hulpmiddel voor het doen van experimenten. De standaardreferentie is het boek

G.H. Golub & C.F. Van Loan, Matrix Computations, The Johns Hopkins University Press, Baltimore, Maryland, USA, 2^dedruk, 1988.

1 een mini-inleiding in MATLAB

“Matlab” is een interactieve taal, ontworpen door Cleve Moler, die is ontstaan uit een demon- stratieproject, waarmee studenten op een eenvoudige wijze zouden kunnen experimenteren met rekenmethoden in de lineaire algebra met gebruik van de standaardroutines uit LINPACK (pakket voor het oplossen van lineaire stelsels vergelijkingen) en EISPACK (pakket voor het eigen- waardeprobleem) zonder zelf de basisalgoritmen gedetailleerd te moeten implementeren. Het ontwerp was zo succesvol, dat Moler een onderneming oprichtte (MathWorks) die het ontwerp heeft uitgebouwd tot een zeer krachtige programmeeromgeving met een uitgebreide verzameling numerieke en grafische hulpmiddelen voor het oplossen (en simuleren) van problemen en het grafisch voorstellen van de oplossing.

De basis-datastructuur is de matrix. Met de instructie “p=5; q=7; A = rand(p,q)” wordt een reële matrix gecreëerd met 5 rijen en 7 kolommen (in IR^5×7 dus) bestaande uit uniform op [0 , 1] verdeelde random getallen. Een matrix met één kolom is een kolomvector, een matrix met

´e´en rij is een rijvector en een enkele “real” staat gelijk met een 1×1–matrix; de typen real en vector zijn dus geen aparte datatypen. Er wordt gerekend met de standaard IEEE 64-bits reals.

Bewerkingen met matrices en vectoren volgen de standaard regels van de lineaire algebra. Een matrix kun je altijd vermenigvuldigen met een scalair (= 1×1–matrix), twee matrices van gelijke afmetingen kun je optellen en aftrekken en als A ∈ IR^p×q en B ∈ IR^q×r, dan is het product de matrix AB ∈ IR^p×r. Het accent wordt gebruikt voor transpositie, dus als A ∈ IR^p×q dan A^′ ∈ IR^q×p. Het systeem rekent even gemakkelijk met complexe matrices volgens de geldende regels uit de lineaire algebra. Een complex getal wordt altijd voorgesteld in de vorm re¨eel deel + i × imaginair deel. Het accent geeft dan transpositie plus complex toevoegen aan.

Voorbeeld(>> is de matlab-prompt):

>> x=[1+i,1-i]

x =

1.0000 + 1.0000i 1.0000 - 1.0000i

>> x’

ans =

1.0000 - 1.0000i 1.0000 + 1.0000i

>> x’*x ans =

2.0000 0 - 2.0000i

0 + 2.0000i 2.0000

>> x*x’

ans = 4

>>

Er zijn allerhande mogelijkheden om deelmatrices te selecteren, b.v. als A ∈ C^p×q, dan is real(A)∈ IR^p×q het re¨ele deel en imag(A)∈ IR^p×q het imaginaire deel. A(:,k) is de k-de kolom (mits 1 ≤ k ≤ q) en A(1:3,2:2:q) is een matrix bestaande uit de elementen met even kolomindex uit de eerste 3 rijen van A. De opdracht x = A\b lost het stelsel vergelijkingen Ax = b op met de optimale methode, dus met Gausseliminatie met rijverwisseling als A vierkant en goed geconditioneerd is, en met een QR-ontbinding of een Singuliere-waardenontbinding als A slecht geconditioneerd of niet vierkant is. De dimensies van b en A moeten natuurlijk compatibel zijn.

Het hele scala van standaard matrix- en vectorroutines, zoals FFT, QR, LU, Choleski, SVD en de berekening van eigenwaarden/eigenvectoren is beschikbaar.

Het systeem is op het eerste gezicht “commandline-oriented” en interpreterend, maar er zijn uitgebreide mogelijkheden tot het oproepen van routines (m-files zonder parametersubstitutie en zonder locale variabelen) en functies (met call by value invoerparameters en uitvoerparameters).

(3)

De standaard controlestructuren zijn (zie help if, help for, etc)

if <boolean> , <statements> else <statements> end for <variable> = <range> , <statements> end

while <boolean> , <statements> end

Een “ , ” (komma) en een “ ; ” (puntkomma) en een <eoln> (end-of-line) worden gebruikt om statements van elkaar te scheiden. Een opdracht in het command-window wordt ge¨ınterpreteerd en onmiddellijk uitgevoerd en het resultaat wordt in het commandwindow afgedrukt, tenzij de statement wordt afgesloten met een “ ; ”.

Omdat iedere programmaregel ge¨ınterpreteerd wordt, is de berekening van het inproduct for k = 1:100 , s = s + a(k) ∗ b(k) ; end

van de kolomvectoren a en b van lengte n, veel trager dan de opdracht s = a’ * b . Ook de statement [m,i]= max(abs(A(1:n,k))) is veel sneller dan het equivalente

m=0; for j = 1:n, if abs(A(j,k)) > m, i = j; m = abs(A(j,k)); end, end . Een declaratie van een matrix of vector is in principe niet nodig. Een datastruktuur wordt automatisch gecre¨eerd en/of uitgebreid, indien nodig. B.v. in

for k = 1:100, c(k) = a(k) ∗ b(k) ; end

wordt de vector c in iedere slag met één element verlengd (en dus worden alle reeds aanwezige elementen in iedere slag gecopieerd naar een nieuwe plaats in het geheugen). Veel sneller zal de volgende variant zijn, waarin vooraf de gehele vector ineens wordt gecreëerd:

c = zeros(100,1); for k = 1:100, c(k) = a (k) ∗ b(k) ; end Nog sneller kunnen we dit doen met een “elementsgewijze” vermenigvuldiging¹ “ . ∗ ”:

c = a .* b

Let wel op, met “c = a .* b” wordt een nieuwe matrix van de gevraagde afmetingen gecre¨eerd, terwijl met “c(:,k) = a .* b” de k-de kolom van de bestaande matrix c wordt overschreven met nieuwe waarde. In dit laatste geval moet deze matrix c wel de goede kolomlengte hebben;

als er reeds k of meer kolommen zijn wordt de k-de overschreven en als er minder zijn, wordt de matrix uitgebreid met het benodigde aantal kolommen (met nullen).

Je kunt in Matlab interactief opdrachten intypen, maar (wat veel handiger is) je kunt ook script-en function- files (in de vorm van m-files) maken, waarin je een “programma” neer- schrijft en die je telkens weer kunt oproepen. Beide files zijn text-files met “extensie” .m . In een script m-file kun je een aantal opdrachten zetten die door matlab worden uitgevoerd als je de naam van de file intypt in het command window; met “echo on” krijg je alle opdrachten achtereenvolgens op je scherm en met “echo off” zet je dit weer af. Alle variabelen die je in een script-file definieert zijn globaal en blijven bestaan na afloop. Het kan dus handig zijn zo’n script file te beginnen met “clear all, clc” om alle aanwezige globale variabelen weg te gooien en het command window leeg te maken, zodat je programma geen gebruik kan maken van rommel, die na een vorige verwerking is achtergebleven.

Een functie (een functie-m-file) is een onafhankelijk programmaonderdeel, dat je kunt oproepen, waaraan je parameters kunt meegeven (invoerparameters) en waarvan je de resultaten via uitvoerparameters weer ter beschikking kunt stellen van het oproepende programmaonderdeel. Voor iedere functie maak je een aparte m-file met naam “<functienaam>.m” en eerste regel “[a,b] = function <functienaam>(x,y,z)”, waarbij “x, y, z” de invoerparameters (call by value) en

“a, b” de uitvoerparameters zijn. Het aantal actuele in- en uitvoerparameters mag kleiner zijn dan het aantal formele. Binnen de functie is het aantal actuele in- en uitvoerparameters op te vragen met “nargin” resp. “nargout”. De nietgespecificeerde invoerparameters kun je dan een

1De vermenigvuldigingsoperator “*”, de delingsoperator “/” en de machtsverheffing “^” volgen de gebruikelijke rekenregels voor matrices. De elementsgewijze operatoren “.*”, “./” en “.^” voeren de bedoelde operatie uit op elementen met gelijke indices, mits beide operanden (matrices) gelijke afmetingen hebben. Als “A” een matrix is, dan is “A*A” of “A^2” het product van “A” met zichzelf en dit product is alleen gedefinieerd als “A” vierkant is; “A.*A”

of “A.^2” is een matrix van gelijke afmetingen als “A” met als elementen de kwadraten van de overeenkomstige elementen van “A”.

(4)

1 EEN MINI-INLEIDING IN MATLAB 4

default-waarde geven. De dimensies van een variabele zijn op te vragen met “size”. Binnen een functie zijn alle variabelen locaal en kan er geen beroep worden gedaan op globaal aanwezige variabelen (behalve “pi” en “i”); dit is dus anders dan bij een oproep van een script-m-file, waarbij de letterlijke inhoud van de file ter plaatse wordt gesubstitueerd en uitgevoerd. Om de naam van een functie als parameter door te geven moet een kunstgreep uitgehaald worden die veel lijkt op het “call by name”-mechanisme van “algol60”. De oproep van de doorgegeven functie wordt als rij karakters opgebouwd en aan de functie “eval” aangeboden ter evaluatie, zoals in het volgende voorbeeld voor het doorgeven van een functienaam aan een functie-m-file (alle karakters op een regel vanaf het %-teken zijn commentaar):

function [y,aantal]=integraal(fun,a,b,n)

% oproep I=integraal(’mijnfun’,0,1,23)

% voor het integreren van de functie met naam "mijnfun".

% Deze functie kan een standaardfunctie zijn in matlab zoals "sin" of "atan"

% of "exp" of een zelfgeschreven m-file met naam "mijnfun.m", zie beneden

% de functie moet zo geschreven zijn dat deze een vector als

% parameter kan meekrijgen en een even grote vector als resultaat aflevert h=(b-a)/n;

x=a+(0:n)*h;

aantal=n+1;

funexpr=[fun,’(x)’]; % opbouw van de functie-oproep z=eval(funexpr); % uitvoering d.m.v. eval

y=(z(1)+z(n+1)+2*sum(z(2:n)))*h/2;

% Voorbeeld van een functiefile:

%% function w=mijnfun(x);

%% w=1./(x.^2+1)

% Als je m-file-functie alleen kan rekenen met een scalaire

% parameter, dan moet je de laatste drie regels van de integraal-m-file

% vervangen door:

%% z=zeros(size(x));

%% for k=1:n+1,funexpr=[fun,’(x(k))’];z(k)=eval(funexpr);end

%% y=(z(1)+z(n+1)+2*sum(z(2:n)))*h/2;

Een alternatieve mogelijkheid in een nieuwere versie van matlab is de mogelijkheid om de functie via @<functienaam> door te geven en de functiewaarden binnen de oproepende functie te evalueren met feval(<functienaam>,<parameters>).

function [y,aantal]=integraal(fun,a,b,n)

% oproep I=integraal(@mijnfun,0,1,23)

% voor het integreren van de functie met naam "mijnfun".

% Deze functie kan een standaardfunctie zijn in matlab zoals "sin" of "atan"

% of "exp" of een zelfgeschreven m-file met naam "mijnfun.m", zie boven

% de functie moet zo geschreven zijn dat deze een vector als

% parameter kan meekrijgen en een even grote vector als resultaat aflevert h=(b-a)/n;

x=a+(0:n)*h;

aantal=n+1;

z=eval(fun,x); % evalueer de functie "fun" met de vector "x" als argument y=(z(1)+z(n+1)+2*sum(z(2:n)))*h/2;

De Matlab-software is voorzien van uitgebreide 2D en 3D grafische hulpmiddelen (met interactieve mogelijkheden, zie “help plot” en “help surf”) en goede documentatie via een “help

<commandnaam>” op de command line en via een help-window. Er is een uitgebreide verzameling

(5)

“toolboxen” voor ijle matrices, signaalverwerking, niet-lineaire vergelijkingen, spline-interpolatie, simulaties, wavelets, gewone en parti¨ele differentiaalvergelijkingen en er is zelfs een interface met Maple.

N.B. 1: “i”, “eps” en “ans” zijn standaard variabelen, de imaginaire eenheid√

−1 , de standaard machineprecisie resp. het laatst berekende en in het command window afgedrukte anonyme antwoord. Je kan problemen verwachten als je deze herdefinieert.

N.B. 2: Onder Linux kun je matlab opstarten door in een terminal window of console na de linux-prompt het commando “matlab&” in te typen. De terminal window kun je daarna min- imaliseren, maar je mag hen niet afsluiten, anders verdwijnt matlab weer. De computers in het computerzaaltje van de vakgroep wiskunde op 10F gebruiken een licentie van een server bij TW;

het opstarten kan nogal eens lang duren.

N.B. 3: De uitkomst van een statement wordt op het scherm afgedrukt, tenzij de statement wordt afgesloten met een “;”. Het formaat van de print wordt bestuurd door de opdracht format met parameters short, long, short e, long e, +, zie help format.

Met de opdracht fprintf heb je volledige controle over de manier waarop de gewenste uitvoer op het scherm (of in een file) komt te staan. Voorbeeld:

x= 0 : .1 : 1; y= [x; exp(x)]; fprintf(^′ %6.2f %20.8d \n^′, y);

Tussen de aanhalingstekens staat het “format” waarmee de inhoud van de vector y wordt afgedrukt, “%6.2f” geeft aan dat het eerstvolgende getal in de invoerrij moet worden afgedrukt in fixed-point notatie op een veld van 6 posities met twee cijfers achter de decimale punt, “%20.8d”

geeft aan dat het eerstvolgende getal in de invoerrij moet worden afgedrukt in floating-point notatie op een veld van 20 posities met een mantisse met acht cijfers achter de decimale punt en

“\n” geeft aan dat de printer moet overgaan naar de volgende regel alvorens volgende getallen af te drukken. De elementen van de vector y worden rijsgewijs afgelopen. Het format-voorschrift wordt telkens herhaald totdat de invoerrij leeg is.

Een uitgebreidere inleiding tot Matlab, A Matlab Primer van Kermit Sigmon, staat op mijn website, http://homepages.vub.ac.be/∼pdegroen/numeriek/matlab primer.pdf .

De zoekwoorden ‘matlab tutorial’ geeft op het internet een grote verzameling van goede inleidin- gen tot Matlab, o.a. die van de producent “Mathworks” zelf.

(6)

2 VOORBEELDEN VAN INSTABIELE ALGORITMEN 6

2 Voorbeelden van instabiele algoritmen

2.a Recursieve berekening van een integraal.

Definieer de integraal

En :=

Z ₁

0

xⁿe^x−1dx voor n = 0, 1, 2, 3, · · · . De waarde van E0 is eenvoudig uit te rekenen,

E₀= Z 1

0 e^x−1dx = e^x−1

1

0 = 1 − e⁻¹ = 0.63212055882856 . Voor andere waarden van n leiden we via parti¨ele integratie de volgende recursie af:

E_n = Z ₁

0

xⁿe^x−1dx = xⁿe^x−1

1 0 − n

Z ₁

0

xⁿ⁻¹e^x−1dx = 1 − nEn−1. De voorwaartse recursie,

E₀ := 1 − e⁻¹ , En:= 1 − nEn−1 (n = 1, 2, · · ·),

is instabiel, zoals blijkt uit de negatieve waarde voor n = 18 in de tabel hieronder. De reden is, dat een fout ε in E_k−1 een fout kε in Ek veroorzaakt. De fout in E18 is dus ongeveer 18! ≈ 10¹⁶ maal die in E₀.

De achterwaartse recursie,

kies E_m= willekeurig , E_n−1 = (1 − En)/n (n = m, m − 1, · · ·),

is stabiel ongeacht de startwaarde voor Em, en geeft de correcte waarde van En als m maar voldoend groot wordt gekozen t.o.v. n . Dit zien we in kolom 4, waar E₁₈ = 0 gekozen is; hier neemt de fout bij iedere iteratie verder af en bij E₅ is deze in de afrondfout verdwenen.

n voorwaarts achterwaarts achterwaarts verschil tussen vanaf k = 0 vanaf n = 50 vanaf n = 18 kolommen 3 en 4 0 0.63212055882856 0.63212055882856 0.63212055882856 0.00000000000000 1 0.36787944117144 0.36787944117144 0.36787944117144 0.00000000000000 2 0.26424111765712 0.26424111765712 0.26424111765712 0.00000000000000 3 0.20727664702865 0.20727664702865 0.20727664702865 -0.00000000000000 4 0.17089341188538 0.17089341188538 0.17089341188538 0.00000000000000 5 0.14553294057308 0.14553294057308 0.14553294057308 -0.00000000000000 6 0.12680235656152 0.12680235656153 0.12680235656152 0.00000000000001 7 0.11238350406936 0.11238350406930 0.11238350406934 -0.00000000000004 8 0.10093196744509 0.10093196744559 0.10093196744528 0.00000000000032 9 0.09161229299417 0.09161229298966 0.09161229299250 -0.00000000000284 10 0.08387707005829 0.08387707010339 0.08387707007499 0.00000000002841 11 0.07735222935878 0.07735222886266 0.07735222917515 -0.00000000031248 12 0.07177324769464 0.07177325364803 0.07177324989825 0.00000000374978 13 0.06694777996972 0.06694770257562 0.06694775132275 -0.00000004874714 14 0.06273108042387 0.06273216394138 0.06273148148148 0.00000068245990 15 0.05903379364190 0.05901754087930 0.05902777777778 -0.00001023689848 16 0.05545930172957 0.05571934593124 0.05555555555556 0.00016379037568 17 0.05719187059731 0.05277111916899 0.05555555555556 -0.00278443638656 18-0.02945367075154 0.05011985495809 0 0.05011985495809

(7)

2.b De berekening van de variantie

De variantie van een stel metingen kan berekend worden met twee Mathematisch equivalente formules. Gegeven n metingen {x1, x₂, · · · , xn} van een grootheid X, dan is zijn gemiddelde g en variantie S² gegeven door

g := 1 n

Xn

k=1

x_k , S_n² := 1 n − 1

Xn

k=1

(x_k− g)² = 1 n − 1

Xn

k=1

x²_k− ng²

! .

De tweede formule is potentieel numeriek instabiel (als S_n² << g²) en veel gevoeliger voor kleine storingen van het gemiddelde g, zoals we zien in het volgende experiment.

Experiment(met Matlab, >> is de matlab-prompt)

>> format short e

>> RelStoringG=1e-12 RelStoringG =

1.0000e-012

>> n=10000;

>> x=randn(n,1)+1e8*ones(n,1);

>> g=sum(x)/n;

>> sig2=x’*x-n*g*g;

>> sig1=(x-g*ones(size(x)))’*(x-g*ones(size(x)));

>> g=sum(x)/n*(1+RelStoringG);

>> sig2s=x’*x-n*g*g;

>> sig1s=(x-g*ones(size(x)))’*(x-g*ones(size(x)));

>> Waarden=[sig1,sig2,sig1s,sig2s]

Waarden =

9.7946e+003 -8.1920e+004 9.7946e+003 -2.0008e+008

>> sprintf([’berekende waarde met formule 1 : %25.15e\ n’,...

’berekende waarde met formule 1 en relatieve storing van g : %25.15e\ n’,...

’berekende waarde met formule 2 : %25.15e\ n’,...

’berekende waarde met formule 2 en relatieve storing van g : %25.15e\ n’]...

,sig1,sig2,sig1s,sig2s) ans =

berekende waarde met formule 1 : 9.794567005712350e+003 berekende waarde met formule 1 en relatieve storing van g : 9.794567105990183e+003 berekende waarde met formule 2 : -8.192000000000000e+004 berekende waarde met formule 2 en relatieve storing van g : -2.000814080000000e+008 De met formule 2 berekende som van kwadraten is in dit experiment zelfs (toevallig) negatief en zeer gevoelig voor een kleine relatieve storing van het berekende gemiddelde.

(8)

3 FOUTENANALYSE 8

3 Foutenanalyse

3.a Elementaire definities

Gegeven is een grootheid X en haar benadering X. De absolute en relatieve fouten in de be-^e naderingX worden gegeven door:^e

absolute fout in X : Fê _X := X − Xê zodat X = X + Fê _X, relatieve fout in X : fê _X := X − Xê

X zodat X = X(1 + f^e _X) (mits X 6= 0).

(3.1)

Het begrip “absolute fout” heeft in principe niets te maken met absolute waarden; absoluut staat slechts in tegenstelling tot relatief. De absolute fout heeft dezelfde dimensies (bv. lengte, gewicht, tijd) als de grootheid zelf, terwijl de relatieve fout dimensieloos is.

Opgave 1: laat zien, dat voor de absolute en relatieve fouten in de twee groothedenX en^e Y geldt:^e F_X+Y = FX+ FY en fX∗Y = fX+ fY + fXfY .

Als we de (absolute of relatieve) fout in een (gemeten of berekende) grootheidX kennen, danê kennen we ook de grootheid zelf! Helaas zijn we bijna nooit in deze situatie en kennen we alleen een bovengrens voor de absolute waarde van de fout. In het gangbare spraakgebruik spreken we gewoonlijk over de “fout” in een grootheid terwijl we zo’n bovengrens bedoelen (of nog erger, terwijl we de spreiding in de stochastische fluctuaties rond de exacte waarde bedoelen). Dus, bij een gegeven benaderingX van een grootheid X definiëren we:ê

∆_X is (een bovengrens voor) de absolute fout inXê als |X − X | ≤ ∆ê X, δ_X is (een bovengrens voor) de relatieve fout inXê als

X − Xe X

≤ δX. (3.2) Opgave 2: Bewijs de volgende rekenregels voor “de fouten” in de grootheden X en^e Y :^e

∆X±Y ≤ ∆^X + ∆Y ∆XY ≤ |Y |∆^X + |X|∆^Y + ∆X∆Y , δ_X±Y ≤ |X|δX + |Y |δY

|X ± Y | δ_XY ≤ δX + δ_Y + δ_Xδ_Y . (3.3) N.B. Lees deze regels alsvolgt: Als ∆X en ∆Y bovengrenzen zijn voor de fouten in X resp. Y , dan is er een bovengrens ∆_X±Y voor de fout in X ±Y waarvoor geldt ∆X±Y ≤ ∆X+∆_Y . Hieruit volgt dus dat ∆X + ∆Y een bovengrens voor de fout in X ± Y is, etc.

Wat zijn de overeenkomstige rekenregels voor de absolute en relatieve fouten (bovengrenzen) in het quoti¨ent X/Y ?

3.b Voorstelling van re¨ele getallen en “floating-point” aritmetiek

Om een groot dynamisch bereik mogelijk te maken voor re¨ele getallen worden deze in een computer opgeslagen in de vorm mantisse maal exponent. Hiertoe wordt een grondtal β (meestal 2, soms 8 (vroeger op CDC) of 16 (IBM)) gekozen. Een x ∈ IR kan dan worden voorgesteld door een paar (m, e) met

x = m · β^e, (3.4)

waarin m de mantisse is en e de exponent. Omdat het paar (m · β , e − 1 ) hetzelfde getal voorstel kunnen we de mantisse normaliseren, b.v. door 1/β ≤ | m | < 1. Het spreekt vanzelf dat we in de praktijk een eindige representatie willen hebben en dus het aantal β-tallige cijfers in

(9)

mantisse en exponent zullen beperken. De IEEE-standaard voor 64-bits REALs is een tweetallige representatie (β = 2) met 53 resp. 10 bits voor de absolute waarden van mantisse en exponent en twee tekenbits. Omdat een genormaliseerde binaire mantisse altijd begint met een 1 (ga na!), hoeft dit eerste bit niet opgeslagen te worden. Met 10 bits is ook de grootte van de exponent aan een maximum gebonden. Getallen die een exponent groter dan 2¹⁰ of kleiner dan 2⁻¹⁰ vragen (waarvan de absolute waarde dus kleiner dan (ongeveer) 10⁻³⁰⁰ of groter dan 10³⁰⁰ is), kunnen dus niet gerepresenteerd worden; we spreken dan van over- of underflow. De IEEE-standaard geeft de mogelijkheid om by underflow een getal op nul te zetten, en bij overflow een N aN (Not a Number) te genereren zodat er een soepele foutenopvang mogelijk is. Een re¨eel getal binnen het bereik zal in het algemeen niet exakt gerepresenteerd kunnen worden. Voor een gegeven x ∈ IR (binnen het bereik) noteren we met fl(x) het meest naburige wel repesenteerbare getal (machinegetal). Het verschil x − fl(x) is dan de afrondfout.

Stelling. Als voor een machinegetal een β-tallige representatie wordt gekozen met t bits in de mantisse, dan geldt voor de relatieve afrondfout bij afronding naar het dichtstbijzijnde machinegetal (behoudens over- en underflow):

x − fl(x) x

≤ η maar ook

x − fl(x) fl(x)

≤ η met η := ¹₂β^1−t. (3.5) De grootheid η heet de machineprecisie.

Opgave 3: Bewijs deze stelling.

Ga ook na, dat er (behoudens over- en underflow) getallen ε1 en ε2 zijn bij iedere aritmetische operatie ⊙ ∈ {+, −, ∗, /} tussen twee machinegetallen x en y, zodat

fl(x ⊙ y) = (x ⊙ y)(1 + ε¹) = x ⊙ y

1 + ε₂ met |ε¹| ≤ η en ε²| ≤ η . (3.6) Opmerking. We kunnen η ook defini¨eren als het grootste re¨ele getal, zodat fl(1 + η) = 1, ga na!

Opgave 4: De reeksontwikkeling van de exponentiaal is: e^x = X∞

k=0

x^k k!

Hoeveel termen heb je nodig om e⁻⁵ te berekenen met een relatieve fout kleiner dan 10⁻³? Kun je dit doen met een computer, waarin de variabelen van het type IR een mantisse van 4 decimalen hebben ?

Is er een betere manier om e⁻⁵ te berekenen met zo’n computer ?

Opgave 5: Laat f een voldoend gladde re¨ele funktie (b.v. f (x) = sin(x)) zijn met maxx |f^′′′(x)| ≤ M .

De afgeleide van f in x kunnen we dan benaderen met de centrale differentie Dhf (x) := f (x + h) − f(x − h))

2h .

Laat zien, dat voor de afbreekfout in D_hf geldt:

f (x + h) − f(x − h)

2h = f^′(x) + h²

6 f^′′′(x + ϑh) met |ϑ| ≤ 1 . (3.7) Veronderstel, dat er voor de berekening van f een procedure beschikbaar is, die bij iedere waarde van x een resultaat aflevert met een relatieve fout kleiner dan of gelijk aan 2η . Geef dan een (goede) bovengrens voor de relatieve fout in de berekende waarde van D_hf als funktie van h en schets een grafiek van (een bovengrens voor) de totale fout (afbreek- plus afrondfout) in deze berekende waarde.

(10)

3 FOUTENANALYSE 10

3.c De onvermijdelijke fout

We beschouwen het probleem om voor een gegeven (gladde, minstens C²) functie f op IR en voor een gegeven argument x de waarde van f (x) te berekenen met een of andere formule of algoritme.

A priori weten we dat we bij berekening op en computer de gegeven x zullen moeten afronden tot een binair “machinegetal” en dat we dus eigenlijk een ander sommetje oplossen, namelijk y = f (x + ξ) met (onbekende) relatieve fout | ξ/x | ≤ η kleiner dan de machineprecisie. Mete Taylor zien we

y = f (x + ξ) = f (x) + ξ fe ^′(x) + O(ξ²) zodat y − y ≈ ξ f_e ^′(x) . Voor de relatieve afwijking geldt dus

y − ye y ≈ ξ

x

x f^′(x)

f (x) en bij benadering y − ye

y

≤ C η met C :=

x f^′(x) f (x)

. (3.8) De foutvermenigvuldigingsfactor C noemen we het conditiegetal van het probleem.

Omdat we het resultaat van de berekeningy nog moeten converteren naar een decimaal antwoord_e y =b y (1+ϑ) met | ϑ | ≤ η vinden we de bovengrens (C +1)η voor de relatieve fout in de berekende_e waarde y . We noemen deze fout de onvermijdelijke fout in de berekende waarde van y._b

3.d Voorbeelden van een afrondfoutenanalyse Gevraagd te berekenen x = ϕ(a).

Met een algoritme voor het berekenen van ϕ(a) vinden we ten gevolge van afrondfouten de berekende waarde: fl(x).

In een foutenanalyse proberen we fouten δx, δa of εa en εx te vinden, zodat fl(x) = x + δ_x voorwaartse foutenanalyse

= ϕ(a + δa) achterwaartse foutenanalyse

= ϕ(a + εa) + εx gemengde foutenanalyse

Definitie: We noemen de algoritme numeriek stabiel als we kunnen bewijzen:

δx of εx van de grootteorde van de onvermijdelijke fout, δa of εa van de grootteorde van de machineprecisie.

Voorbeeld 1: Er is een ε met | ε | ≤ η ( = machineprecisie ) zodat

fl(a + b) =





a + b + ε (a + b) voorwaarts

ea +^eb met _ea := a(1 + ε) en ^eb := b(1 + ε) achterwaarts Voorbeeld 2: Er zijn ε₁ en ε₂ ( met | εi| ≤ η ) zodat

fl(1 − x²) = (1 − x ∗ x ∗ (1 + ε¹)) ∗ (1 + ε²)

= (1 − xe²) (1 + ε₂) met x := x_e √

1 + ε₁ gemengd.

Voorbeeld 3: Geef een schatting van de afrondfout in de berekende waarde van de positieve wortel van de vierkantsvergelijking

a − 2x − c x² met a ≥ 0 en c ≥ 0 bij gebruik van de formule

x := −1 + √

1 + a c c

(11)

onder de aanname betreffende de afrondfout in de berekende waarde van de vierkantswortel fl(√

x) = √

x (1 + ε_x) met | εx| ≤ η voor iedere x . Antwoord: Er bestaan ε₁, ε₂ en ε₃ met | εⁱ| ≤ η , zodat

fl(√

1 + a c) = ^p(1 + a c (1 + ε₁)) (1 + ε₂) (1 + ε₃)

= √

1 + _eac (1 + ξ₁) met ξ₁ := √

1 + ε₂(1 + ε₃) − 1 en a := a (1 + εe ₁) Bijgevolg zijn er ξ₂ en ξ₃, ( | ξⁱ| ≤ η ) zodat:

fl(x) = −1 + √

1 + _ea c (1 + ξ1)

c (1 + ξ2) (1 + ξ3)

= −1 + √

1 + _ea c

c (1 + ξ₂) (1 + ξ₃) +

√1 + _ea c

c ξ₁(1 + ξ₂) (1 + ξ₃) De tweede term kan groot zijn t.o.v. x als | ac | ≪ 1.

Alternatieve (numeriek stabiele) rekenwijze voor deze wortel:

x := a

1 + √

1 + a c.

Voorbeeld 4: Afrondfout in de berekende waarde van het inprodukt S :=

Xn

i=1

x_iy_i berekend met algoritme: S := 0;

fori := 1 to n do S := S + x_i ∗ yi

Voor de berekende waarde van S vinden we getallen ξi en ηi met | ξⁱ|, | εⁱ| ≤ η, i = 1 · · · n : fl(S) = x₁y₁(1 + ξ₁) (1 + ε₂) · · · (1 + εn)

+ x2y2(1 + ξ2) (1 + ε2) · · · (1 + εⁿ) + · · ·

+ x_n−2y_n−2(1 + ξ_n−2) (1 + ε_n−2) · · · (1 + εⁿ) + x_n−1y_n−1(1 + ξ_n−1) (1 + ε_n−1) (1 + ε_n) + xnyn(1 + ξn) (1 + εn)

zodat

S − fl(S) = Xn

i=1

x_iy_iζ_i met

ζi := 1 − (1 + ξⁱ) (1 + εi) · · · (1 + εⁿ) en | ζⁱ| ≤ (n − i + 2) η als nη ≤ 0.1 . Bijgevolg geldt voor de voorwaartse fout:

|S − fl(S)

S | ≤ (n + 1) η

| S | Xn

i=1

| xiy_i| ≤ (n + 1) ηk x k2k y k2

| x^Ty| (3.9)

Voorbeeld 5: Bereken xnuit de vergelijking a =

Xn

i=1

xiyi, a, x1 · · · xⁿ⁻¹, y1 · · · yⁿ gegeven,

(12)

3 FOUTENANALYSE 12

en bepaal de afrondfout in de berekende waarde van xn. Algoritme:





S := a;

fori := 1 to n − 1 DO S := S − xⁱ ∗ yⁱ; xn := S / yn

Voor de berekende waarden van S en xn vinden we voor zekere ξi en ηi met | ξⁱ|, | εⁱ| ≤ η : fl(S) = a (1 + ε₁) · · · (1 + εn−1)

− x1y₁(1 + ξ₁) (1 + ε₁) · · · (1 + εn−1)

− x2y₂(1 + ξ₂) (1 + ε₂) · · · (1 + εn−1)

− · · ·

− xn−2y_n−2(1 + ξ_n−2) (1 + ε_n−2) (1 + ε_n−1)

− xn−1y_n−1(1 + ξ_n−1) (1 + ε_n−1) en

xe_n := fl(x_n) = fl(S) / ( y_n(1 + ξ_n) )

Deling door (1 + ε₁) · · · (1 + εn−1) geeft de achterwaartse foutschatting:

a = x1y1(1 + ξ1) + x2y2

1 + ξ₂

1 + ε1 + · · · + x_n−1y_n−1 1 + ξ_n−1

(1 + ε₁) · · · (1 + εn−2) + x_enyn

1 + ξ_n

(1 + ε₁) · · · (1 + εn−1)

=

n−1X

i=1

xiyi(1 + δi) + x_enyn(1 + δn) met

δ_i := 1 + ξ_i

(1 + ε₁) · · · (1 + εi−1) − 1 , zodat | δi| ≤ (i + 1) η als n η < 0.1 . Conclusie: De berekende waarde x_e_n is de oplossing van een naburige vergelijking

a = Xn j=1

xjyej, y_ej := yj(1 + δj) . (3.10)

Voorbeeld 6, foutschatting voor En. In hoofdstuk 2a van de syllabus hebben we afgeleid:

E_n= 1 − n En−1. (3.11)

LaatE^en:= fl(En) de berekende waarde van Enzijn, dan zijn er getallen ξnen ζnwaarvoor geldt:

Ee_n = fl(1 − fl(nE^e_n−1)) = (1 − nE^e_n−1(1 + ξ_n))/(1 + ζ_n) , | ξn| ≤ η en | ζn| ≤ η , (3.12) of anders:

Een+ ζnEen= 1 − nE^e_n−1− n ξⁿE^e_n−1. (3.13) Trekken we hiervan vergelijking (3.11) af, dan vinden we een recursie voor de verschillen met de exacte waarden

Een− Eⁿ= −n(Eê_n−1− En−1) − ζⁿEên− n ξⁿEê_n−1 (3.14) of met Fn:=Eên− Eⁿ en δn:= −ζⁿEên− n ξⁿEên−1

Fn= −n Fn−1+ δn, F₀ = fl(E₀) − E0, | F0| ≤ η E0 ≤ η. (3.15)

(13)

Merk op dat uit de positiviteit van En en vergelijking (3.11) volgt, dat En−1 ≤ 1/n en dat dit ook voorE^e_n−1 moet gelden zolang dat een redelijke benadering van Enis; onder deze voorwaarde geldt dus ook | δⁿ| ≤ 2η . Fⁿ voldoet dan dus aan de ongelijkheid

| Fⁿ| ≤ n | Fn−1| + 2η . (3.16)

Bijgevolg is er een majorerende rij {F^bn} zodat

| Fn| ≤F^b_n met F^b_n= nF^b_n−1+ 2η , F^b₀= η . (3.17) De recursie voorF^bn levert een a priori bovengrens voor de fout in de berekende waarde van En:

|E^en− Eⁿ| = Fⁿ ≤F^bn= n! η

1 + 2

1!+ 2

2!+ · · · + 2 n!

≤ n! η (2 e¹− 1) . (3.18) Bij n = 18 is deze bovengrens reeds veel groter dan 1, zodat dan waarschijnlijk al niet meer voldaan is aan de conditie |E^e_n| ≤ 1.

Een betere bovengrens kunnen we berekenen, door tegelijk metE^en ook een bovengrens voor de fout in deze berekende waarde uit te rekenen. Uit (3.14) volgt immers

| Fn| ≤ n | Fn−1| + η |Eê_n| + n η |Eê_n−1| (3.19) Als we dus tegelijk metEê_n de meelopende foutschatting Fê_n berekenen met de recursie

Fe_n= nFê_n+ η |Eê_n| + n η |Eê_n−1| (3.20) dan vinden we zo voor iedere n a posteriori een (vrij goede) bovengrens voor de absolute fout in de berekende waarde van Enuit algoritme (3.11). Merk op dat we deze (kleinere) bovengrens pas na de berekeningen kunnen kennen, omdat deze beter rekening houdt met feitelijk gemaakte afrondfouten.

3.e oefeningen

1. Vorm de volgende expressies om tot een numeriek stabiele vorm 1

1 + 2x−1 − x

1 + x voor | x | ≪ 1 (3.21)

r x + 1

x − r

x − 1

x voor | x | ≫ 1 (3.22)

1 − cos x

x voor | x | ≪ 1 (3.23)

2. We zouden de functie x 7→ arctan(x) kunnen evalueren via de identiteit arctan x = arcsin x

√1 + x² (3.24)

Bepaal (een goede bovengrens voor) de relatieve afrondfout in het resultaat, als we rekenen met een machineprecisie η, en bepaal voor welke waarden van x de genoemde methode acceptabel is.

3. Laat f een C³-functie zijn op IR waarvan de derde afgeleide begrensd is door 1: | f⁽³⁾(x) | ≤ 1 . We berekenen de afgeleide met de centrale differentie

f^′(a) = f (a + h) − f(a − h)

2h + O(h²) . (3.25)

Geef een schatting van de totale fout in het resultaat (afrond– + afbreekfout) als we rekenen op een computer met machineprecisie η en als we een routine voor de berekening van f (x) ter beschikking hebben die de waarde ervan aflevert met een absolute fout van ten hoogste ε.

(14)

4 LINEAIRE ALGEBRA 14

4. Een polynoom P van de graad n kan voorgesteld worden door een som of een product P (x) :=

Xn

k=0

akx^n−k of P (x) := a₀ Yn

k=1

(x − x^k) (met a₀ 6= 0 ),

met co¨effici¨enten a0, a1, · · · , aⁿ resp. (complexe) nulpunten x1, x2, · · · , xⁿ en met a₀ 6= 0 . De waarde van het polynoom wordt voor de gegeven waarde ξ van x het best berekend met de algoritme van Horner:

b₀ := a₀; for k := 1 to n do b_k := b_k−1∗ ξ + ak end; dan volgt: P (ξ) = b_n . (3.26) Desgewenst kan de waarde D van de afgeleide P^′(ξ) tegelijkertijd meeberekend worden,

D := 0 ; P := a₀; for k := 1 to n do D := D ∗ ξ + P ; P := P ∗ ξ + a^k end. a. Bewijs correctheid van de algoritme (3.26).

b. Laat zien, dat b₀, · · · , bn−1 de co¨effici¨enten zijn van het polynoom P (x)/(x − ξ) van de graad n−1 als ξ een nulpunt van P is; dit wil zeggen, dat we met de algoritme van Horner een nulpunt kunnen uitdelen (synthetische deling).

c. Laat zien dat er getallen δkzijn zodat de berekende waarde fl(P (x)) gelijk is aan de exacte waarde van een naburig polynoom,

fl(P (x)) = Xn

k=0

ea_kx^n−k met _ea_k:= a_k(1 + δ_k) en | δk| ≤ (2n − 2k + 1)η + O(η²) .

d. Laat zien, dat de volgende algoritme een “meelopende fout” d berekent,

P := a⁰; d := 0; for k := 1 to n do d := d + | P | ; P := P ∗x+a^k ; d := d ∗ | x | + | P | end zodat na afloop van de algoritme geldt:

| fl(P (x)) − P (x) | ≤ d η .

5. Voor de standaardafwijking S bestaan in de statistiek twee formules die wiskundig (in exacte re¨ele arithmetiek) gelijkwaardig zijn :

S² = 1 n − 1 (

Xn

i=1

x²_i − n g² ) en S² = 1 n − 1

Xn

i=1

(xi− g)² met g het gemiddelde:

g := 1 n

Xn

i=1

xi .

Welke van de twee zou je gebruiken in een numeriek programma en waarom?

4 Lineaire Algebra

4.a notaties

Hoewel het voor bewijzen in de de lineaire algebra vaak handig is om met een abstracte vectorruimte E van dimensie n over het lichaam van de re¨ele of complexe getallen te werken, is het rekenen met vectoren en afbeeldingen alleen mogelijk op basis van een representatie van vectoren en afbeeldingen ten opzichte van een basis in de ruimte. Dit wil zeggen dat een vector in IRⁿ of Cⁿ altijd een kolommetje van n re¨ele of complexe getallen is en dat een lineaire afbeelding

(15)

van een n-dimensionale ruimte naar een m-dimensionale altijd wordt gegeven door een matrix in IR^m×n of C^m×n.

• Een vector x ∈ IRⁿ is een kolommetje van n re¨ele getallen,

x=





 x₁ x₂ ... x_n





 met componenten x1, · · · , xⁿ. (4.1)

In druk noteren we x met een vetgedrukte letter en in manuscript met een ondersteepte letter x .

• Een matrix A ∈ IR^m×n noteren we altijd met een hoofdletter en de matrixelementen aij met een kleine letter. De kolommen van een matrix zijn vectoren in IR^m; zij spannen tesamen de beeldruimte Im(A) op:

A =





a11 · · · a1n

... . .. ... am1 · · · a^mn



= (a₁| a2| · · · | an) , zodat ak=



 a_1k

... amk



. (4.2)

De overeenkomstige notatie in matlab is: als A een matrix is, dan is de vector A(:,k) de k–de kolom ervan.

Een vector kunnen we opvatten als een matrix met ´e´en kolom; ook matlab doet dat.

• Een matrix A ∈ IR^m×n en een vector x ∈ IRⁿ kunnen we alsvolgt partitioneren:

A =

A₁₁ A₁₂ A₂₁ A₂₂

en x =

x₁ x₂

zodat A x =

A₁₁x₁+ A₁₂x₂ A₂₁x₁+ A₂₂x₂

(4.3) mits de dimensies kloppen:

A₁₁∈ IR^p×r, A₁₂∈ IR^q×r, A₂₁∈ IR^p×s, A₂₂∈ IR^q×s u∈ IR^p, v ∈ IR^q, p + q = n en r + s = m

De overeenkomstige notatie in matlab is: als A een matrix is, dan wordt het deel A₂₂weergegeven door B=A(p+1:n , r+1:m) . Denk er wel aan dat de indices in B dan verschoven zijn, zodat B(1,1)=A(p+1,r+1)etc.

• De getransponeerde van een matrix A wordt gegeven door A^T; als de matrix complex is, dan wordt de Hermitisch geconjungeerde (transponeren en complex toevoegen) gegeven door A^H. In matlab wordt A^H verkregen door een accent: A’ .

• De norm van een vector x ∈ IRⁿ wordt gegeven door k x k . Zoals bekend zijn alle normen in een eindigdimensionale vectorruimte equivalent (wat is dat? waarom?). In het vervolg worden vrijwel uitsluitend de Euclidische norm k · k2 (of ℓ²-norm) , de max-norm k · k∞ (of ℓ^∞-norm) en de 1–norm k · k1 (ℓ¹–norm of duale van de max-norm) gebruikt:

k x k1:=

Xn j=1

| x^j| , k x k²2 :=

Xn j=1

| x^j|² en k x k∞:= max

j | x^j| (4.4) Bij de Euclidische norm hoort een inproduct:

Als u , v ∈ Cⁿ, dan hu , vi := u^Hv= Xn

j=1

ujvj. (4.5)

Aangezien u^H een rijvector of een 1×n–matrix is, kunnen we het inproduct (4.5) ook aanduiden met het matrix–matrix product u^Hv(of u^Tvvoor re¨ele vectoren) en met u’*v in matlab.

De vectoren u , v heten orthogonaal als u^T v= 0 .

(16)

• De matrixnorm A 7→ k A k , geassoci¨eerd met de vectornorm x 7→ k x k , wordt voor een matrix A ∈ IR^m×ngedefini¨eerd door

k A k := max

x∈ IRⁿ, kxk 6= 0

k A x k

kxk = max

x∈ IRⁿ, kxk = 1 k A x k (4.6)

In de teller staat er natuurlijk de vectornorm in IR^m en in de noemer die in IRⁿ. De zo gedefini¨eerde matrixnorm wordt ook wel een lub–norm genoemd (lub van least upper bound).

Ga na dat een lub–norm altijd voldoet aan de producteigenschap k A B k ≤ k A k k B k .

Voor de lub–normen bij de normen (4.4) gebruiken we dezelfe notatie k · k1, k · k2 en k · k∞.

• De Frobenius–norm van een matrix A ∈ IR^m×n is gedefinieerd door k A k²F :=

Xm

i=1

Xn

j=1

| a^ij|². (4.7)

Deze norm voldoet wel aan de producteigenschap k A B kF ≤ k A kFk B kF maar is geen lub–norm.

In feite is het de Euclidische norm van de matrix gezien als element van een m×n–dimensionale vectorruimte.

In matlab kun je de genoemde vector– en matrixnormen voor een object a eenvoudig berekenen met de functie norm(a,p) waar p een van de symbolen “1”, “2”, “inf” of “’fro’” is (de laatste met quotes!).

• Een (re¨ele) matrix A ∈ IR^n×n heet orthogonaal als A^TA = I de identitieke afbeelding in IRⁿ is;

een (complexe) matrix A ∈ C^n×nheet unitair als A^HA = I. Gan na, dat dit impliceert A A^T = I resp. A A^H = I.

Als A ∈ IR^m×n met m > n en A^TA = I , dan heeft A orthonormale kolommen en dan noemen we A een parti¨ele isometrie.

• Een diagonaalmatrix D ∈ IR^m×nis een matrix waarvan alle elementen buiten de hoofddiagonaal nul zijn, dus D = (dij) met dij = 0 als i 6= j .

Voor een gegeven vector a ∈ IRⁿ defini¨eren we de diagonaalmatrix D := diag(a) ∈ IR^m×n met m ≥ n door dii= a_i; we kiezen m = n, tenzij uit de context duidelijk is, wat de waarde van m is.

De matlab-functie diag construeert uit een vector een vierkante matrix met deze vector als hoofddiagonaal; Als de functie op een m×n–matrix wordt toegepast, extraheert deze de diagonaal en levert een vector af met lengte min(m,n).

4.b oefeningen

1. Bewijs voor de geassocieerde matrixnormen de volgende identiteiten:

kAk1 = max

j

Xn

i=1

|a^ij| , kAk∞= max

i

Xn

j=1

|a^ij| en kAk2 = max

x6=0,y6=0

|(Ax, y)|

kxk2kyk2

met (x, y) :=^Pⁿ_i=1x_iy_i.

2. Bewijs voor x ∈ IRⁿ en A ∈ IR^n×n de ongelijkheden:

1) kxk2 ≤ kxk1≤√

nkxk2 2) kxk∞≤ kxk2 ≤√ nkxk∞

3) 1

√nkAk2 ≤ kAk1 ≤√

nkAk2 4) 1

√nkAk∞≤ kAk2 ≤√

nkAk∞

Toon met voorbeelden aan dat de ongelijkheden scherp zijn, d.w.z. dat er gelijkheid kan optreden.

3. Laat zien dat de 2-norm van een matrix unitair invariant is (kU Ak2 = kAk2 voor iedere unitaire afbeelding U ).

(17)

4. Toon aan dat de “Frobenius”-norm niet geassocieerd is aan een vectornorm. Bewijs er de producteigenschap voor en laat ook zien dat deze norm invariant onder unitaire transformaties (kU Ak^F = kAk^F voor iedere unitaire afbeelding U ).

5. Toon voor de Frobeniusnorm aan:

k A k²F = spoor ( A^TA ) = som der eigenwaarden van A^T A k A k²2 = grootste eigenwaarde van A^TA

√1

nk A k^F ≤ k A k² ≤ k A k^F

N.B. De wortels van de eigenwaarden van A^TA zijn de singuliere waarden van A.

6. Voor een gegeven vector a ∈ IRⁿ is de afbeelding fa:= x 7→ a^Tx een lineaire transformatie van IRⁿnaar IR. Bewijs dat k fak1 = k a k∞, k fak∞= k a k1 en k fak2 = k a k2.

4.c De singuliere-waardenontbinding

Theorem: Bij iedere (re¨ele) matrix A ∈ IR^m×n zijn er orthogonale matrices U ∈ IR^m×m en V ∈ IR^n×n en p := min{m, n} niet-negatieve getallen σ1, · · · , σ^p zodat

A = U Σ V^T, Σ := diag(σ₁, · · · , σp) ∈ IR^m×n. (4.8) NB. De getallen σ₁, · · · , σ^p heten de singuliere waarden van A.

Het is gebruikelijk om de singuliere waarden te ordenen in een dalende rij, zodat σ_k≥ σk+1. In matlab kan de singuliere-waardenontbinding berekend worden met de functie svd:

s=svd(A)levert de singuliere waarden af in de vector s.

[U,S,V]=svd(A)levert in U, S en V de drie matrices van de ontbinding (4.8) af.

Bewijs. We kunnen dit op twee manieren bewijzen, met en zonder de eigenwaardenontbinding van A^TA . We kiezen voor de eenvoud m ≥ n. We gebruiken natuurlijk steeds de Euclidische vectornorm en de geassocieerde matrixnorm.

1. De matrix A^TA ∈ IR^n×n is symmetrisch en niet-negatief definiet. Zij heeft dus n niet- negatieve eigenwaarden, die we dalend kunnen ordenen λ₁≥ λ2 ≥ · · · ≥ λⁿ ≥ 0 . Verder is er een bijbehorende orthonormale basis van eigenvectoren v₁, · · · , vⁿ zodat A^TAv_k= λ_kv_k. Definieer nu

Σ := diag(√

λ₁, · · · ,√

λ_n) , V := (v₁| · · · | vn) en U :=^e

Av₁

√λ₁ · · ·

Avn

√λn

,

dan is V een orthogonale matrix en heeft de matrixU orthonormale kolommen. We kunnen deze^e matrix aanvullen met m−n kolommen tot een orthogonale matrix (hoe?). Het resultaat voldoet dan aan (4.8).

2. Een elegant elementair bewijs gaat met volledige inductie alsvolgt. Definieer σ₁ := kAk . De functie x 7→ kAxk is continu en neemt op de eenheidsbol {k x k = 1} zijn maximum aan.

Er is dus een vector v₁ met norm 1 zodat kAv1k = kAk = σ1. Definieer u₁ := Av₁/σ₁ en maak orthogonale matrices U := (u₁|U ) en V := (v^b ₁|V ) met deze vectoren als eerste kolom^b door de stelsels {u¹} en {v¹} aan te vullen tot orthonormale bases in IR^m resp. IRⁿ (b.v. met Gram-Schmidt). Hiermee vinden we:

V := (v₁|V )^b zodat A V = A (v₁|V ) = (Av^b ₁| AV ) = (σ^b ₁u₁| AV )^b en

A := Ue ^TA V = u^T₁ Ub

!

(σ₁u₁| AV ) =^b σ₁u^T₁u₁ u^T₁AV^b σ₁U^b^Tu₁ U^b^TAV^b

!

= σ₁ w^T 0 A^b

!

. (4.9)

(18)

Hierbij is 0 de nulvector, omdat de kolommen vanU bij definitie loodrecht staan op u^b ₁.

De (rij-)vector w^T := u^T₁AV is de eerste rij van^b A behalve het eerste element; we zullen bewijzen^e dat ook dit een nulvector is. De overblijvende matrixA :=^b U^b^TAV ∈ IR^b (m−1)×(n−1)is van kleinere dimensie. We kunnen nu uitrekenen:

Ae σ₁

w

!

2

=

σ₁²+ w^Tw Awb

!

2

≥ (σ1²+ w^Tw)²,

maar, omdat de matrixnorm invariant is onder orthogonale transformaties geldt ook

Ae σ₁

w

!

2

≤ σ²1

σ₁ w

!

2

= σ₁²(σ₁²+ w^Tw) .

Bijgevolg geldt σ₁²+ w^Tw≤ σ1² en moet de vector w ∈ IRⁿ⁻¹ de nulvector zijn. Het gevolg is, dat

U^TA V = σ₁ 0^T 0 A^b

!

. (4.10)

We kunnen hetzelfde argument dus toepassen op de kleinere matrixA en zo doorgaan totdat deze^b triviaal is.

De singuliere-waardenontbinding (Engels: singular value decomposition, afgekort SVD) is een zeer geschikt hulpmiddel om de rang van een matrix ( = dimensie van de beeldruimte) te bepalen, je telt gewoon het aantal singuliere waarden, dat niet nul is. In de praktijk stelt zich echter het probleem, dat de verzameling inverteerbare matrices (of, als m 6= n, matrices van volle rang) dicht liggen in de verzameling van alle matrices, zodat willekeurig dicht bij een matrix van lagere rang nog matrices van volle rang liggen. Ten gevolge van afrondfouten bij numerieke berekeningen kun je dus nooit beslissen of een matrix defectief is (niet van volle rang). Daarom defini¨eren we de numerieke rang rang(A, ε) als de rang van de matrix met kleinste rang in een bol met straal ε rond A:

rang(A, ε) := min {rang(A + E) | E ∈ IR^m×n, kEk ≤ ε} (4.11) Als we de 2-norm gebruiken en de singuliere waarden σ₁, · · · , σ^p van A berekenen met p = min{m, n}, dan vinden we

rang(A, ε) := r als σ₁ ≥ · · · ≥ σ^r > ε ≥ σr+1≥ · · · ≥ σ^p, (4.12) want, als A = U Σ V^T, dan voldoet E := U diag(0 , · · · , 0 , σr+1, · · · , σp) V^T aan de eis kEk ≤ ε en rang(A − E) = r .

Voorbeeld van data reductie met de SVD. Het magische vierkant uit de houtsnede “Me- lencolia” van D¨urer is gescand in een 359×371 matrix van grijswaarden. Het linker plaatje van fig. 2 toont een grafiek van de singuliere waarden van deze matrix. Het rechter plaatje toont de grijswaardenverdeling die we overhouden als we in de SVD alle singuliere waarden op nul zetten behalve de grootste. Kennelijk is het rooster het dominante element in dit plaatje. Het middelste plaatje krijgen we als we alleen de 36 grootste singuliere waarden behouden en alle andere (de kleintjes) op nul zetten. Kennelijk levert dit al een zeer goede benadering van het originele plaatje.

Deze analyse techniek is niet erg gebruikelijk voor plaatjes. In de statistiek is dit echter een zeer gebruikelijke vorm van datareductie, genaamd “principal component analysis”. Het rooster in het rechter plaatje van fig. 2 is de eerste “hoofdas” van de grijswaardenmatrix van het plaatje van het magisch vierkant; Het linker plaatje is de “scree² plot”.

2Scree is het naar beneden gerolde puin dat onder aan een zeer steile helling ligt.

(19)

Duerer, Melancholia 359 x 371 pixels

Figure 1: D¨urer’s houtsnede “melencolia” (links) en als detail

“het magische vierkant” in de rechter bovenhoek (rechts).

0 50 100 150 200 250 300 350

10⁰ 10¹ 10² 10³ 10⁴

singular values, logarithmic plot reduction to q = 36 singular values reduction to q = 1 singular values

Figure 2: Grafiek van de 359 singuliere waarden van de grijswaardenmatrix van het detail (links) en reconstructies ervan met 36 (midden) en precies 1 singuliere waarden (rechts).

4.d Het Conditiegetal

Voor een gegeven matrix A ∈ IR^n×n en vector b ∈ IRⁿ willen we methoden ontwikkelen om het stelsel lineaire vergelijkingen

A x = b (4.13)

op te lossen. Alvorens hier aan te beginnen is het nuttig om de gevoeligheid van het probleem voor kleine verstoringen van A en b te bestuderen. Naast het “exacte” probleem (4.13) beschouwen we dus het verstoorde probleem

(A + E) (x + w) = b + d , met E ∈ IR^n×n en d ∈ IRⁿ klein. (4.14) en we proberen een schatting te vinden voor de resulterende afwijking w in de oplossing. Hiervoor hebben we de volgende schatting nodig voor de norm van (A + E)⁻¹:

Lemma. Als A ∈ IR^n×ninverteerbaar is en als voor E ∈ IR^n×ngeldt kA⁻¹Ek < 1 , dan is A + E inverteerbaar en voldoet aan de schatting

k(A + E)⁻¹k ≤ kA⁻¹k

1 − kA⁻¹Ek (4.15)

Bewijs. Als I de identiteit in IRⁿ is en F ∈ IR^n×n voldoet aan kF k < 1 , dan geldt voor alle x∈ IRⁿ dat

kIx + F xk ≥ kxk − kF xk ≥ (1 − kF k)kxk > 0 .

Geen enkele niet-triviale vector wordt dus op de nulvector afgebeeld, zodat de kern van I + F de nulruimte is en I + F inverteerbaar is. Voor alle x ∈ IRⁿ geldt dus

kxk = k(I + F )(I + F )⁻¹xk ≥ (1 − kF k)k(I + F )⁻¹xk zodat k(I + F )⁻¹k ≤ 1 1 − kF k.

Numerieke Lineaire Algebra

Contents

1 een mini-inleiding in MATLAB

2 Voorbeelden van instabiele algoritmen

3 Foutenanalyse

4 Lineaire Algebra