Minstakvadratmetoden i matlab

Antag att vi har en matris A som har fler rader än kolumner och en kolumnvektor b som är lika hög. Då är Ax=b ett överbestämt system.

Kort härledning av minstakvadratmetoden

Vi har det överbestämda systemet Ax=b. Låt oss studera kolumnvektorerna i A och kalla dem för A1, A2, ... An. På samma sätt kallar vi de olika komponenterna i x för x1, x2, ... xn. Då gäller att A1*x1 + A2*x2 + ... + An*xn = b. Försäkra dig gärna om att detta stämmer!

Med andra ord vi försöker uttrycka b med hjälp av en linjärkombination av kolumnvektorerna i A. Vi har dock fler ekvationer än obekanta, systemet är överbestämt.

Låt oss säga att vi har m ekvationer och n obektanta. Detta betyder att vi har ett b som befinner sig i ett m-dimensionellt rum, men vi har bara n stycken basvektorer (kolumnvektorerna i A) för att försöka uttrycka b med. Problem!

Det är oftast trevligt när man kan tolka problem geometriskt. Här gäller att kolumnvektorena i A spänner upp ett hyperplan och vi ska försöka hitta den punkten på hyperplanet som ligger så nära b som möjligt.

Jag upprepar, vi försöker uttrycka b som en linjärkombination av kolumnvektorerna i A, detta kan dock inte göras exakt:

Ax = A1*x1 + A2*x2 + ... + An*xn = b

Vad är felet vi gör när vi väljer en punkt Ax på hyperplanet?

Felet är residualen, dvs r = Ax - b. När är felet som minst? Jo - när residualen är ortogonal mot hyperplanet, för då befinner vi oss i den punkten på hyperplanet som är så nära den sökta punkten b som möjligt. När är residualen ortogonal mot hyperplanet? Det är den när den är ortogonal mot alla kolumnvektorerna i A, ty dessa spänner upp hyperplanet.

På mattespråk: A'*r = A'*(Ax-b) = 0.

Flyttar vi om i ekvationen får vi den bekanta formen på minstakvadratmetoden: A'Ax = A'b.

Problem som kan uppstå, illa konditionerat

Överbestämda system kan lösas genom att projicera kolumnerna i både A och b på A:s kolumnrum och sedan lösa det nya systemet som fås. Problemet är dock att om A har ett högt konditioneringstal så kommer A'*A ha A:s konditioneringstal i kvadrat!

Det vill säga, att be matlab lösa (A'*A)\(A'*b) kan leda till följande:

Detta går dock att kringå genom QR-faktorisering, men hur gör man det? Den goda nyheten är att matlab gör det automagiskt om ni skriver A\b.

Vad är då QR-faktorisering och hur hjälper det oss?

När man QR-faktoriserar en matris så skriver man om den på formen Q*R där Q och R är två matriser. Det påminner lite om när man skriver om ett komplext tal på polär form. Tänk er Q som vinkeln och R som beloppet ("längden"). Q har normen ett och Q*R=A, vidare så är Q'*Q=I.

Det ursprungliga problemet A*x=b kan skrivas Q*R*x=b där vi bara bytt ut A mot QR.

Vi vill minimera residualen e, vilket är felet vi gär då vi försöker uttrycka b med bara A:s kolumnvektorer.

|e| = |A*x-b| = |Q*R*x-b| = /Q har normen ett, ändrar ej värdet / = |Q'*Q*R*x - Q'*b| = /vi vet att Q'*Q = I/ = |R*x-Q'*b|

Dvs e är som minst då R*x = Q'*b, här tar vi inte A'*A så vi slipper problemet med stort konditioneringstal. Det är egentligen såhär matlab gör när den löser minstakvadratproblem.

De goda nyheterna som sagt är att det enda vi behöver göra för att matlab ska lösa minstakvadratproblemet på det smarta sättet är att skriva:

Tillbaka till numbio04 (email)