5 Totale kleinste kwadraten

5 TOTALE KLEINSTE KWADRATEN

5

5.a

49

Totale kleinste kwadraten

Beste benadering in IR

Als we de verzameling punten V := {x1 , x2 , · · · , xm } in IR hebben gegeven en we vragen welk punt z het dichtst bij al deze punten ligt, dan is het antwoord natuurlijk afhankelijk van de definitie van “het dichtstbij”. Het eenvoudigste criterium is het minimaal zijn van de “de som van kwadraten van de verschillen”; we moeten dan een z zoeken zodat J(z) =

m X

(xi − z)2

i=1

(5.1)

minimaal is. Deze functionaal is een nette parabool in z met als uniek minimum het gemiddelde P x := m i=1 xi , zoals we zien uit de ongelijkheid J(z) =

m X i=1

2

(xi −z) =

m X i=1

2

2

{(xi −x) +(x−z) +2(xi −x)(x−z)} ≥

m X i=1

(xi −x)2 = J(x) ∀ z ∈ IR ,

omdat de som van de dubbele produkten nul is. We zien dus, dat de kleinste-kwadratenbenadering van de verzameling punten V in IR gegeven wordt door het gemiddelde x .

5.b

Lineaire regressie in IR2

Een grootheid y is afhankelijk van een variabele x en we vermoeden een lineair verband, y = a + bx,

(5.2)

zoals bijvoorbeeld de stroom door een weerstand bij een gegeven spanningsverschil of de uitrekking van een veer als funktie van het gewicht dat we er aan hangen. We doen een aantal metingen van y als funktie van x, resulterend in de koppels {(xi , yi )|i = 1, · · · , m}. In het ideale geval liggen deze koppels dus op de lijn gegeven door (5.2) voor zekere waarden van a en b. Ten gevolge van meetfouten en andere afwijkingen van het ideaal, liggen de meetpunten echter niet precies op een rechte en stelt zich het probleem om uit de gegeven metingen {(xi , yi ) | i = 1, · · · , m} parameters a en b te bepalen, zodat de lijn y = ax + b zo goed mogelijk er bij past, zie fig. 9a. Als we x als de onafhankelijke variabele beschouwen (die bij een meting dus vooraf gekozen kan worden) dan hebben we per meting een fout yi − a − bxi . We gaan dan a en b zo kiezen dat de “totale” fout minimaal is. Een gebruikelijke keuze voor deze totale fout is de som van de kwadraten (die voor de analyse verreweg het eenvoudigst is) J(z) :=

m X i=1

2

(yi − a − bxi )

met

z :=

Ã !

a ∈ IR2 b

(5.3)

¡ ¢

en het probleem is dan een koppel z := ab ∈ IR2 te bepalen, dat J(z) minimaliseert. In formule (5.3) lijkt het probleem tweedimensionaal. Als we echter in IRm de vektoren x, y en e en de matrix A definiëren,    

x := 

x1 x2 .. . xm



  , 

   

y := 

y1 y2 .. . ym



  , 

 

1 1   e :=  ..  . 1

en



1 1  A :=  .. .

1

x1 x2 .. . xm



   ∈ IRm×2 , 

(5.4)


5

5

x x

0

50

x

x

x

x

0

x

x x

x x

x

x

x

x x

x x

x

-5

x

x

x

-5 0

0

a. beste rechte met y als funktie van x

b. beste rechte met x als funktie van y

5

0

5

x x x

x

x

0

x

x x

x x

x

x

x

x

-5

x x

x x

x

x

x

x

-5

0 c. totale kleinste-kwadratenbenadering

0 d. de drie benaderende rechten tesamen

Figure 9: Datapunten met de beste bijpassende rechten; (a) geeft de beste rechte als de som van kwadraten van de vertikale afstanden (gestippeld) wordt geminimaliseerd, (b) als de som van kwadraten van horizontale afstanden wordt geminimaliseerd en (c) als de som vam kwadraten van de afstanden van de punten tot de lijn wordt geminimaliseerd. In (d) zijn de drie “benaderingen” tesamen afgebeeld. dan is som van kwadraten J(z) = lengte van de verschilvektor

Pm

i=1 (yi

− a − bxi )2 precies het kwadraat van de (Euclidische) zodat J(z) = ky − Azk22 .

y − ae − bx = y − Az

De lengte van de verschilvektor y − Az is minimaal (zie fig. 10) als deze loodrecht staat op de beeldruimte Im(A) van A; d.w.z.: y − Az ⊥ vect{e, x} ofwel, met (·, ·) als notatie voor het inprodukt, µ

¶

w , y − Az = 0

∀ w ∈ Im(A).

Aangezien Im(A) = {A w | w ∈ IR2 } geldt dus µ

¶

Aw , y − Az = 0

∀ w ∈ IR2 .

We mogen A getransponeerd naar de andere zijde van het inprodukt overbrengen, zodat µ

T

T

¶

w , A y − A Az = 0

∀ w ∈ IR2


51

Az

©©

©©

©©

Im(A) ©©

© *© © © © © A ©© A ©© © © Ay − Az ©© © © © A © © © © A © © © © -AU © ©

O

y

Figure 10: De vector y, zijn (orthogonale) projektie op Im(A) en de residuvektor y − Az . en zo vinden we de normaalvergelijkingen voor het minimalisatieprobleem (5.3): AT Az = AT y .

(5.5)

Dit heten de normaalvergelijkingen omdat het residu y − Az normaal is (= loodrecht staat) op Im(A). De verkregen rechte y = a + bx wordt (in de statistiek) de “regressierechte” van y op x genoemd. Een eenvoudiger weg naar het minimum van (5.3) vinden we, als we de data verschuiven zodat hun gemiddelde nul is. In plaats van (5.2) werken we dan met het (equivalente) model y − y = α + β(x − x),

y :=

m 1 X yi , m i=1

x :=

m 1 X xi . m i=1

(5.6)

In dit geval moeten we de funktie e J(α, β) :=

n X i=1

(yi − y − α − β(xi − x))2 = ky − ye − αe − β(x − xe)k22

minimaliseren. Omdat de vektor e loodrecht staat op de vektoren y − ye en x − xe is Je minimaal als Pm (xi − x)(yi − y) (x − xe)T (y − ye) Pm = i=1 . (5.7) α = 0 en β = T 2 (x − xe) (x − xe) i=1 (xi − x) Hieruit zien we, dat de regressierechte (5.2) of (5.6) door het zwaartepunt (x, y)T van de puntenwolk gaat.

5.c

Lineaire regressie van x op y

Analoog aan het model (5.2) of (5.6), waarin y als funktie van x is beschouwd, kunnen we ook andersom x als funktie van y beschouwen, x − x = γ + δ(y − y) .

(5.8)

Minimalisatie van de som van kwadraten van horizontale afstanden geeft analoog aan (5.7) γ=0

en

δ=

Pn

i=1 (xi − x)(yi − Pm 2 i=1 (yi − y)

y)

.

(5.9)

Omdat we hier de som van kwadraten van horizontale afstanden minimaliseren, zie fig. 9b, vinden we (i.h.a.) een andere regressielijn dan in het voorgaande geval. Het gemeenschappelijke punt is precies het zwaartepunt (x, y)T van de data.


5.d

52

De totale kleinste-kwadratenbenadering

Als er bij de gegevens {(xi , yi )|i = 1, · · · , m}

geen duidelijke voorkeur is voor de keuze van y als funktie van x of van x als funktie van y, d.w.z. als beide grootheden (random) meetfouten bevatten, is het beter om een benaderingscriterium te kiezen, dat geen voorkeur voor een van beide bevat. Het ligt dan voor de hand om som van de kwadraten van de (echte) afstanden van de datapunten tot de lijn te minimaliseren, zie fig. 9c. Uit de analytische meetkunde weten we, dat de verzameling ℓ := {z ∈ IR2 | rT (z − w) = 0} voor gegeven r ∈ IR2 een rechte beschrijft door het punt w loodrecht op de vektor r en dat de afstand van een willekeurig punt x ∈ IR2 tot deze lijn gegeven wordt door de lengte van de projektie van de verschilvektor x − w op een veelvoud van r, zie fig. 11. ℓ = {z ∈ IR2 | (z − w, r) = 0} Z Z

Z

Z

x

o ¶ 7S Z ¶ T S Z ¶dist(x, ℓ) = |r √(x−w)| Z S rT r Z¶ S } Z Z Z S Z Z Z Z Sx−w Z Z S Z Z Z Z S Z Z S ¶ 7r Z Z S ¶ Z Z ¶ Z ZS Z Z Sw ¶ Z ¶ Z ¶ Z

O Figure 11: Getekend zijn de vektor r, de lijn door w loodrecht op r, het punt x, de verschilvektor x − w en de ontbinding ervan langs de lijn en loodrecht erop. Als krk2 = 1, dan wordt de afstand van x tot de lijn ℓ gegeven door |rT (x − w)| . Voor het bepalen van de totale kleinste-kwadratenbenadering moeten we dus de loodvektor r van de rechte (van lengte 1) en een vektor w op die rechte bepalen zodat de functionaal I, I(r, w) :=

m ³ X

´2

T

r (xi − w)

=

m ³ X

´2

r1 (xi − w1 ) + r2 (yi − w2 )

(

met

xi = (xi , yi )T

r = (r1 , r2 )T (5.10) minimaal is. Natuurlijk geeft iedere vektor w op de rechte dezelfde rechte; we zullen eerst laten zien, dat het zwaartepunt van de data op de rechte ligt, zodat de keuze w = (x, y)T altijd goed is en er daarna alleen nog de minimalisatie van (5.10) naar r = (r1 , r2 )T overblijft. Splitsen van een term als (xi − w1 )2 in (5.10) geeft: i=1

m X i=1

(xi − w1 )2 =

m ³ X i=1

i=1

´

(xi − x)2 + (x − w1 )2 + 2(xi − x)(x − w1 )

(5.11)

en hierin is de som van de dubbele produkten nul; het analoge geldt voor de term (yi −w2 )2 . Hieruit volgt: I(r, w) = =

Pm

i=1

Pm

i=1

³

³

´2

r1 (xi − w1 ) + r2 (yi − w2 ) ´2

r1 (xi − x) + r2 (yi − y)

≥ I(r, (x, y)T ) ∀w ∈ IR2 .

³

´2

+ m r1 (x − w1 ) + mr2 (y − w2 )

(5.12)


53

Hieruit volgt, dat bij iedere r de keuze w = (x, y)T de functionaal I(r, · ) minimaliseert, zodat het zwaartepunt op de rechte moet liggen. Om de minimale r te vinden, herschrijven we I alsvolgt,

I(r, (x, y)T ) = kBrk22 = rT B T B r ,

met

   

B := 



x1 − x x2 − x .. .

y1 − y y2 − y .. .

xm − x

ym − y

  . 

(5.13)

Het minimum wordt dus aangenomen, als r de rechter singuliere vektor van B is, die behoort bij de kleinste singuliere waarde (of als r de eigenvektor is, behorende bij de kleinste eigenwaarde van B T B). Dit minimum is uniek, als beide singuliere waarden van elkaar verschillen. Als r2 6= 0 volgt hieruit, r1 (5.14) y = y + ρ(x − x) met ρ := − . r2 We merken op, dat het inzicht, dat het zwaartepunt van de data op de gezochte rechte ligt, de sleutel is in de reduktie van (5.10) naar de bepaling van een singuliere waarde en van de oplossing (5.14).

5.e

Regressie in meer dan twee dimensies

We beschouwen opnieuw m koppels metingen {(xi , yi )|i = 1, · · · , m}. In plaats van het polynoom (5.2) van graad 1 willen we nu een polynoom van graad n−1 (met n > 0) vinden, dat hierbij zo goed mogelijk past in “kleinste-kwadratenzin”. Als het polynoom gegeven is door p(x) = c0 + c1 x + c2 x2 + c3 x3 + · · · + cn−1 xn−1

(5.15)

met onbekende coëfficiënten c0 , · · · , cn−1 , dan kunnen we het benaderingsprobleem formuleren als het zoeken naar het minimum van de functionaal J(c) :=

m X i=1

(yi − c0 − c1 x − c2 x2 − · · · − cn−1 xn−1 )2 = ky − Ack22 ,

(5.16)

waar de vektoren y ∈ IRm , c ∈ IRn en de matrix A ∈ IRm×n gedefinieerd zijn door    

y := 



y1 y2 .. . ym

  , 

   

c := 



c0 c1 .. . cn−1

  , 

en



1 1  A :=  .  ..

x1 x2 .. .

1 xm



· · · xn−1 1  · · · xn−1 2  . ..  . 

(5.17)

· · · xn−1 m

Een gelijkaardig probleem krijgen we, als we voor een (te meten) grootheid y, die van verscheidene parameters afhangt, een n-dimensionaal lineair model postuleren, y = c0 + c1 x1 + c2 x2 + · · · + cn−1 xn−1 .

(5.18)

Deze vergelijking beschrijft een hypervlak van dimensie n−1 in een n-dimensionale ruimte IRn . (i) (i) Uit m (m > n) metingen {(x1 , · · · , xn−1 , yi ) | i = 1, · · · , m} willen we weer de beste coëfficiënten in kleinste-kwadratenzin bepalen door het minimaliseren van de functionaal J(c) :=

m X i=1

met



  y :=  

y1 y2 .. . ym



  , 

(i)

(i)

(yi − c0 − c1 x1 − · · · − cn−1 xn−1 )2 = ky − Ack22 

  c :=  

c0 c1 .. . cn−1



  , 

en



1

 1 A :=  .  ..

(1)

x1

(2)

x1 .. .

(m)

1 x1

(1)

· · · xn−1 (2)

(5.19)  

· · · xn−1  . ..  .  (m)

· · · xn−1

(5.20)


54

We merken op, dat we in dit tweede geval analoog aan (5.6) de data kunnen verschuiven naar het zwaartepunt (x1 , · · · , xn−1 , y) met het equivalente model y − y = c0 + c1 (x1 − x1 ) + c2 (x2 − x2 ) + · · · + cn−1 (xn−1 − xn−1 )

met y :=

(5.21)

m m 1 X 1 X (i) yi en xk := xk . We moeten dan de functionaal Je minimaliseren, m i=1 m i=1

e J(c) :=

m µ X i=1

(i)

¶2

(i)

yi − y − c0 − c1 (x1 − x1 ) − · · · − cn−1 (xn−1 − xn−1 )

e 2 . (5.22) = ky − ye − Ack 2

Aangezien de eerste kolom van Ae gelijk is aan de vektor e, staat deze loodrecht op y − ye en op e Dus is de co¨ alle andere kolommen van A. efficiënt c0 automatisch nul. We zien hieruit, dat het zwaartepunt van de data in het hypervlak ligt, opgespannen door vergelijking (5.18).

5.f

De normaalvergelijkingen in IRm

Laat A ∈ IRm×n met m ≥ n een matrix zijn en b ∈ IRm een vektor. In het algemeen zal b geen element van Im(A) zijn, zodat het probleem Ax = b geen oplossing heeft. Wel kunnen we ons analoog aan (5.2) en (5.3) afvragen, welke x de best bijpassende is in kleinste-kwadraten zin en de funktionaal J(x) := kAx − bk22 (5.23)

minimaliseert. We zoeken dus een y = Ax ∈ Im(A) (in de beeldruimte van A), die een minimale afstand heeft tot b.

Zoals we in fig. 10 zien, wordt dit punt Ax gegeven door de orthogonale projektie van b op Im(A) en dus door de eis b − Ax ⊥ Im(A). Evenals in §5.b leidt dit tot de normaalvergelijkingen AT Ax = AT b

(5.24)

Een alternatief bewijs van de bewering “de vektor x minimaliseert de functionaal J in (5.23) dan en slechts dan als x oplossing is van het stelsel normaalvergelijkingen (5.24) AT Ax = AT b” is alsvolgt: Als J(x) minimaal is, dan moet gelden J(x + w) ≥ J(x) ∀w ∈ IRn . Uitwerking van de norm geeft: kAx + Aw − bk2 = kAx − bk2 + 2(Ax − b, Aw) + kAwk2 ≥ kAx − bk2 .

Hieruit volgt, dat de lineaire term nul moet zijn voor alle w: (Ax − b, Aw) = 0

∀w ∈ IRn

en dit is equivalent met de normaalvergelijkingen (5.24). Dit stelsel normaalvergelijkingen (5.24) kunnen we oplossen met een Choleski-ontbinding van de (symmetrische) matrix AT A, op voorwaarde dat deze matrix inverteerbaar is. Deze matrix is (als zij al inverteerbaar is) echter vaak slecht geconditioneerd tengevolge van de kwadratering van A. Voor de nauwkeurigheid van de numerieke berekeningen is het dan beter om een andere methode te gebruiken, waarbij het niet nodig is om het produkt AT A te berekenen. Dit zijn de zogenaamde orthogonale methoden, die gebruik maken van een QR-ontbinding van A, zie het boek van Golub & Van Loan. Voorbeeld. Als we rekenen op een processor met floating point machineprecisie η (met correcte √ afronding) en als ε = 13 η een klein getal is, zodat voor de berekende waarde van 1 + 2ε2 geldt f l(1 + 2ε2 ) = 1, dan vinden we voor de matrix A ∈ IR3×2 

1  A := 0 0



1 ε  ε

de berekende waarde

T

f l(A A) =

µ

1 1

1 1

¶

,

(5.25)


55

zodat de rang van de p berekende AT A gelijk is aan 1, terwijl de rang van A gelijk is aan 2 en het conditiegetal κ2 (A) ≈ 2/η groot is maar toch nog ver verwijderd is van 1/η . Kennelijk verdwijnt alle informatie over de kleine singuliere waarde volledig bij de berekening van AT A .

5.g

Oplossing via de singuliere-waardenontbinding

Omdat het rechterlid AT b van (5.24) een element is van Im(A), heeft dit stelsel vergelijkingen altijd een oplossing. Als echter de matrix A niet van volle rang is, dan is de oplossing niet uniek. Voor de uniciteit kunnen we dan als extra eis stellen, dat de oplossing van het kleinstekwadratenprobleem (5.23) een minimale lengte (Euclidische norm) moet hebben. Om die oplossing te vinden gebruiken we de singuliere-waardenontbinding (voortaan af te korten met SVD, Singular Value Decomposition) van A, A = U ΣV T , (5.26) met U ∈ IRm×m en V ∈ IRn×n orthogonaal en

Σ = diag(σ1 , σ2 , · · · , σn ) ∈ IRm×n ,

waar de singuliere waarden dalend geordend zijn, σk ≥ σk+1 , en waar σr 6= 0 en σk = 0 voor k > r, zodat de rang van de matrix A gelijk is aan r (r ≤ n). Als we de vektoren uk ∈ IRm en vk ∈ IRn definiëren als de k-de kolommen van U resp. V , zodat µ

U = u1 | · · · | um

¶

en

µ

¶

V = v1 | · · · | v n ,

(5.27)

dan vormen de verzamelingen {u1 , · · · , um } en {v1 , · · · , vn } orthonormale bases in IRm resp. IRn en we kunnen A dan schrijven als A=

r X

σk uk vkT .

(5.28)

k=1

De vektoren x en b kunnen we schrijven als lineaire combinaties van deze basisvektoren, x=

n X

ξk vk met ξk onbekend

en

b=

m X

βk uk met βk := uTk b ,

k=1

k=1

en invullen in formule (5.23); we vinden dan J(x) = k

r X

k=1

σk ξk uk −

m X

k=1

βk uk k22 =

r X

k=1

(σk ξk − βk )2 +

m X

βk2 .

(5.29)

k=r+1

Deze kwadratische expressie is minimaal als ξi = βi /σi voor i = 1, · · · , r ongeacht de waarden van ξi voor i = r + 1, · · · , n. Het is duidelijk dat kxk minimaal is als alle coëfficiënten in deze laatste groep nul zijn. De oplossing van het minimalisatieprobleem (5.23) met kleinste lengte wordt dus gegeven door: r X βk T x= uk b vk . (5.30) σ k k=1

De afbeelding van IRm naar IRn , die b afbeeldt op deze minimum-norm-oplossing van (5.23) heet de pseudoinverse (of Moore-Penrose-inverse) van A en wordt aangeduid met A† . Uit (5.30) vinden we analoog aan (5.29): r X 1 A† = (5.31) vk uTk = V Σ† U T . σ k k=1

De matrix Σ† vinden we uit Σ door transpositie plus het vervangen van de positieve diagonaalelementen σk door 1/σk (k = 1 · · · r).

Formule (5.29) beschrijft precies de minimum-norm-oplossing van (5.23) in termen van de SVD; over de feitelijke berekening van de SVD zullen we het later hebben.


5.h

56

Totale kleinste kwadraten in IRn

In het n-dimensionale regressieprobleem (5.18) kunnen we analoog aan (5.10) ook naar het minimum zoeken van de som van kwadraten van Euclidische afstanden van de datapunten tot het hypervlak. De aanpak van §5.d laat zich hierbij onmiddellijk generaliseren. Het n−1-dimensionale hypervlak in IRn door w loodrecht op r is de verzameling {z | ∈ IRn | rT (z − w) = 0} de afstand van de vektor x tot dit hypervlak is rT (x − w) als krk = 1 . We moeten dus weer vektoren r (van lengte 1) en w bepalen, zodat de functionaal I,

I(r, w) :=

m X i=1

2

T

(r (xi − w)) =

m X i=1

 

n X

(i)

2

rj (xj − wj )

j=1

met



(i)



x  1.   .   . 

xi :=  P

. (i)   xn−1 

(5.32)

yi

m 1 minimaal is. Wegens (5.11) geldt ook hier I(r, w) ≥ I(r, x) , waar x := m i−1 xi het zwaartepunt van de data is, zodat we de vektor r vinden door minimalisatie van I(r, x),

I(r, x) = kBrk22 ,



(1)

x1 − x1

 (2)  x1 − x1

B :=  

met



.. . (m) x1 − x1

(1)

· · · xn−1 − xn−1 (2)

· · · xn−1 − xn−1 .. . (m)

· · · xn−1 − xn−1

y1 − y y2 − y .. . ym − y



  .  

(5.33)

Het minimum wordt dus weer aangenomen door de rechter singuliere vektor behorende bij de kleinste singuliere waarde van B en dit minimum is uniek, als deze kleinste singuliere waarde strikt kleiner is dan de andere. Het kleinste-kwadratenprobleem (5.23) is iets algemener dan de lineaire regressie in (5.19), omdat A niet noodzakelijk een kolom (1, 1, · · · , 1)T bevat zoals in (5.20). Het idee voor de aanpak van totale kleinste kwadraten voor een algemeen overbepaald stelsel Ax = b is echter analoog aan het bovenstaande. In §5.f hebben we de functionaal J, gedefinieerd in (5.23), ge¨ınterpreteerd als afstand in IRm en opgelost door een projektie op Im(A), de deelruimte opgespannen door de kolommen van A. Analoog aan 5.d kunnen we Ax = b ook interpreteren in IRn+1 als het zoeken van een n-dimensionaal hypervlak, dat het best past bij een wolk van m punten (observaties). Hiertoe beschouwen we de m rijen van de uitgebreide matrix ( A | −b ) ∈ IRm×(n+1) als vektoren in IRn+1 , 

ak1 ak2 .. .

   ak :=    akn

−bk

      

zodat

( A | −b )T = ( a1 | a2 | · · · | am ) .

(5.34)

De k-de rij ak van ( A | −b ) bevat dus precies de coördinaten van het k-de punt van de puntenwolk (de k-de observatie). De functionaal J kunnen we nu interpreteren als de som

J(x) =

m X

k=1

b )2 (aTk x

met





x1  µ ¶   x2  x  ..  b := =  .  ∈ IRn+1 . x   1  xn  1

(5.35)

b de k-de component is van Ax − b, sommeert J precies de kwadraten van het Aangezien aTk x overschot in iedere component.


57

b gelijk is aan 1, kunnen we de grootheid aTk x b ook interOmdat de n+1-ste component van x b⊥, preteren als de afstand gemeten langs de n+1-ste co¨ ordinaatas van het punt ak ∈ IRn+1 tot x b . Dit betekent, dat het “foutenmodel” d.w.z. tot het hypervlak door de oorsprong loodrecht op x J ervan uit gaat, dat alle fouten in de benadering aan de n+1-ste coördinaat (het rechterlid dus) moeten worden toegeschreven. In de praktijk bevatten meestal alle componenten van een waarneming en dus ook de elementen van A meetfouten en is het natuurlijker de totale fout in alle richtingen te verdisconteren.

De totale kleinste-kwadratenbenadering voor het oplossen van het overbepaalde stelsel Ax = b b ⊥ te bestaat er dan ook in, de Euclidische afstanden van de observaties ak (k = 1 · · · m) tot x ⊥ T b wordt gegeven door ak x b /kx b k. We moeten dus een vektor minimaliseren. De afstand van ak tot x x zoeken, die de functionaal I minimaliseert, m X b k2 b )2 k( A | −b ) x (aTk x = I(x) := bT x bT x b b x x

met

k=1

b := x

µ

x 1

¶

.

(5.36)

b k2 / x bT x b kunnen beperken tot de (comAangezien we ons bij de minimalisatie van k( A | − b ) x b k = 1, is er altijd een minimum x b ; dit geeft echter alleen een minimum voor I als pacte) bol kx b niet gelijk is aan 0. Precies als in (5.33) is I minimaal, als x b de de laatste component van x rechter singuliere vektor is van de matrix ( A | − b ) behorende bij de kleinste singuliere waarde σmin . Als σmin strikt kleiner is dan alle andere singuliere waarden, dan is de oplossing uniek. Als de laatste component van de singuliere vektor gelijk is aan nul, is er geen oplossing voor het totale kleinste-kwadratenprobleem. Er is dus duidelijk een verschil tussen het n-dimensionale lineaire-regressieprobleem (5.18), waar er altijd een oplossing bestaat (in de zin van totale kleinste kwadraten), en het algemene kleinste-kwadratenprobleem (5.23), waar dit niet noodzakelijk het geval is.

Voorbeeld 1. Zoek de regressierechte door de punten (1, 1), (−1, 1), (1, −1) en (−1, −1) , de vier hoekpunten van een vierkant in het platte vlak. Het zwaartepunt is de oorsprong; alle regressierechten gaan door dit punt. Regressie van y op x. Uit formule (5.7) volgt, dat de richtingscoëfficiënt β = 0 in de rechte voor regressie van y op x: de rechte is {(x, y) ∈ IR2 | y = 0 ∀x}.

Regressie van x op y. Analoog volgt uit (5.9), dat de richtingscoëfficiënt δ = 0 in de rechte voor regressie van x op y: de rechte is {(x, y) ∈ IR2 | x = 0 ∀y}. Deze staat dus loodrecht op de vorige. Regressie met totale kleinste kwadraten. Volgens (5.13) moeten we hier de volgende SVD bepalen: 

1  1 B :=   −1 −1



1 −1  = 1  −1



1  2  1  2   1  −2  − 12

1 2 − 12 1 2 − 21

q

1 2

0

0

q

1 2

0



 q  2 1   0 2  q   0 1  2 0

0



0 µ 1 2   0 0  0

0 1

¶

.

(5.37)

Beide singuliere waarden zijn gelijk, zodat iedere rechte door de oorsprong de functionaal I in (5.13) minimaliseert. De totale kleinste-kwadratenbenadering is dus niet uniek; de beide regressierechten zijn het wel (in dit geval). Voorbeeld 2. Los op in (gewone) kleinste-kwadratenzin en in totale kleinste-kwadratenzin: 

1  0 0

  Ã !  µ 1 0 x 1    met normaalvergelijkingen 1 0 = 0 y 1 0

0 0

¶ Ã !

x y

=

µ

1 0

¶

. (5.38)

Uit de normaalvergelijkingen vinden we x = 1 en y is onbepaald, zodat de oplossing voor het (gewone) kleinste-kwadratenprobleem niet uniek is; iedere vektor (1, y)T is een oplossing.


58

Om de totale kleinste-kwadratenbenadering te vinden moeten we van de volgende matrix B de kleinste singuliere waarde bepalen: 

1 B= 0 0

0 0 0



1 1 . 1

Aangezien B rang 2 heeft (B heeft twee onafhankelijke kolommen) en de vektor u := (0, 1, 0)T in de kern van B zit (Bu = 0), is σmin := 0 de kleinste singuliere waarde en u de bijbehorende rechter singuliere vektor. De laatste component van deze vektor is echter 0, zodat herschaling hiervan nooit een 1 kan maken. We zien hier dus een voorbeeld, waarin de functionaal I uit (5.36) een minimum heeft, dat niet overeenkomt met een oplossing van het totale kleinste-kwadratenprobleem.

5.i

Een alternatieve benadering voor totale kleinste kwadraten in IRn

Bij het oplossen van het overbepaalde (en dus i.h.a. strijdige) stelsel vergelijkingen Ax = b in de zin van (gewone) kleinste kwadraten zoeken we een vektor r ∈ IRm van minimale lengte met de eigenschap b + r ∈ Im(A), zodat het gewijzigde stelsel Ax = b + r een oplossing heeft. Alle onzekerheid in de data wordt hierbij toegeschreven aan het rechterlid b. In een benadering via “totale kleinste kwadraten” willen we ook een stuk van de onzekerheid toeschrijven aan de matrix A. We zoeken dus een stoormatrix E ∈ IRm×n en een vektor r ∈ IRm van minimale grootte, zodat b + r ∈ Im(A + E), of anders gezegd, zodat (A + E)x = b + r een oplossing heeft. De grootte van de storing (E | r) ∈ IRm×(n+1) kunnen we meten in de “Frobenius-norm” k·kF (wortel van de som van kwadraten van alle matrixelementen), die evenals de k · k2 -norm invariant is onder orthogonale transformaties. Het probleem kunnen we dus formuleren als het bepalen van het minimum van { k (E | r) kF | E ∈ IRm×n , r ∈ IRm , b + r ∈ Im(A + E) } .

(5.39)

Laten we nu aannemen, dat rang(A) = n, zodat alle kolommen van A onafhankelijk zijn en geen der singuliere waarden van A nul is. Verder nemen we aan, dat b 6∈ Im(A), omdat er anders een oplossing is (met E = 0 en r = 0). Bijgevolg is de rang van (A | b) gelijk aan n + 1. Nu moet (A + E)x = b + r een oplossing hebben. Deze vergelijking kunnen we ook schrijven als b=0 (A + E | b + r) x

met

b := x

µ

x −1

¶

∈ IRn+1 .

Dit betekent, dat de matrix (A + E | b + r) singulier moet zijn en dus, dat de rang ervan hoogstens n mag zijn. Het minimaliserende koppel (E | r) (in de Frobeniusnorm) kunnen we vinden door de SVD van (A | b) te maken, (A | b) = U Σ V T =

n+1 X k=1

σk uk vkT ,

met

 m×m ,   U = (u1 | · · · | um ) ∈ IR

V = (v | · · · | v

) ∈ IR(n+1)×(n+1) ,

1 n+1   Σ = diag(σ , · · · , σ m×(n+1) , 1 n+1 ) ∈ IR

waar U en V orthogonale matrices zijn en Σ een diagonaalmatrix is met de dalend geordende singuliere waarden op de hoofddiagonaal. Het minimaliserende koppel (E0 | r0 ) is dan gelijk aan T (E0 | r0 ) = −σn+1 un+1 vn+1 .

Dit koppel is uniek, als de ongelijkheid σn+1 < σn strikt is. Aangezien (A + E0 | b + r0 ) =

n X

σk uk vkT ,

k=1

zit de rechter singuliere vektor vn+1 , behorend bij de kleinste singuliere waarde, in de kern ervan. b is hier dus een veelvoud van; dit geeft weer de eerder gevonden oplossing van het totale De vektor x kleinste-kwadratenprobleem, als de laatste component van deze vektor niet gelijk is aan nul.

REFERENCES

112

References [1] M. Hestenes & E. Stiefel, Methods of conjugate gradients for solving linear systems, J. Research NBS, 49, pp. 409 – 436, 1952. [2] C. Lanczos, An iteration method for the solution of the eigenvalue problem of linear differential and integral operators, J. Research NBS, 45, pp. 255 – 282, 1950. [3] J.K. Reid, On the method of conjugate gradients for the solution of large sparse systems of linear equations, Proc. Conf. on Large Sparse Sets of Linear Equations, Academic Press, New York, 1971. [4] J.A. Meijerink and H.A. van der Vorst, An iterative solution method for linear systems of which the coefficient matrix is a symmetric M-matrix, Math.of Comp., 31, pp. 148 – 162, 1977. [5] G.H. Golub & C.F. Van Loan, Matrix Computations, The Johns Hopkins University Press, Baltimore, Maryland, USA, 1ste druk, 1983, 2de druk, 1988, 3de druk, 1995. [6] R. Bulirsch & J. Stoer, Introduction to Numerical Analysis, Springer Verlag, Berlin, 1977. (Ook verkrijgbaar in een goedkope duitstalige pocketeditie). [7] D. Kincaid & W. Cheney, Numerical Analysis, Brooks & Cole Publishing Company, Pacific Grove, California, USA, 1991; 2de druk, 1996.

5 Totale kleinste kwadraten

Recommend Documents