Minstakvadratmetoden

Fil:Linear regression.svg

En rät linje som en approximativ beskrivning av en datamängd

Fil:Ickelinjär approximation.png

Icke-linjär approximation (röd kurva) av en punktmängd (blå punkter). Den numeriska modellen är <math>\log(a + b\, x^2)</math> eller <math>\log(1.50632 + 1.42632 x^2)</math> efter tillämpning av minstakvadratmetoden

Minstakvadratmetoden (även minsta-kvadrat-metoden eller minsta kvadrat-metoden) används bland annat vid regressionsanalys för att minimera felet i en funktion som ska anpassas utifrån observerade värden. Exempel på tillämpningar är

• Utifrån gjorda folkräkningar vill man förutsäga befolkningsökningen i ett område genom att göra folkmängden till en funktion av tiden.
• Inom hydrologi vill man beräkna hur stort skyfall som inträffar en gång var hundrade år, till exempel för att kunna dimensionera en mindre damm (se även frekvensanalys). I detta fall görs regnmängden till en funktion av återkomsttiden.

Minstakvadratmetoden har en linjär och en icke-linjär variant beroende på om residualerna (”felen”) är linjära eller inte med avseende på alla obekanta. Den linjära varianten tillämpas inom regressionsanalys och har en sluten form. Den icke-linjära bygger vanligen på iterativa metoder. Vid varje iteration approximeras lösningen med en linjär lösning, varför de grundläggande beräkningarna är snarlika i båda fallen.

Historik

Fil:Carl Friedrich Gauss.jpg

Carl Friedrich Gauss

På nyårsdagen 1801 upptäckte den italienske astronomen Giuseppe Piazzi dvärgplaneten Ceres. Under 40 dagar kunde han följa dess väg, tills Ceres försvann bakom solen. Under året hade många forskare utan framgång försökt att beräkna banan baserat på Piazzis iakttagelser - under antagandet att banan var cirkulär, eftersom endast sådana bandelar kunde bestämmas matematiskt utifrån de observerade positionerna vid denna tidpunkt. Den 24-årige Carl Friedrich Gauss kunde dock beräkna elliptiska banor utifrån tre olika observationer. Med tillgång till betydligt fler spårpunkter, använde han sin minstakvadratmetod för att öka noggrannheten. När småplaneterna på nytt observerades av Franz Xaver von Zach och Heinrich Wilhelm Olbers i december 1801, i exakt de positioner som förutsagts av Gauss, var detta inte bara en stor framgång för Gauss metod, utan ledde även till ett återupprättande av Piazzis rykte, som skadats på grund av konflikten med omloppsbanor beräknade under antagandet att banorna var cirkulära. Minstakvadratmetoden blev snabbt standardförfarandet vid behandlingen av astronomiska och geodetiska mätresultat.

Minstakvadratmetoden tillskrivs vanligen Carl Friedrich Gauss (1795), men publicerades först av Adrien-Marie Legendre.

Anpassning av en funktion till observerade data

En vanlig modell för att representera en mätserie

<math>(x_1,\ y_1),\ (x_2,\ y_2),\,\dots,\ (x_n,\ y_n)</math>

i form av en funktion, är en linjärkombination av m kända (valda) funktioner

<math>f(t) = c_1 f_1(t) + c_2 f_2(t)+\dots + c_m f_m(t)</math>

där koefficienterna c₁, c₂, ... , c_m skall bestämmas för att i minstakvadratmetodens mening bäst anpassa kurvan f till mätserien, vilket innebär att summan

<math>\sum_{i=1}^n\ [y_i-f(x_i)]^2</math>

skall minimeras.

För en lösning konstrueras först den så kallade designmatrisen

<math>A = \begin{bmatrix}

f_1(x_1) & f_2(x_1) &\cdots &f_m(x_1) \\ f_1(x_2) & f_2(x_2) &\cdots &f_m(x_2) \\ \vdots &\vdots &\vdots &\vdots \\ f_1(x_n) & f_2(x_n) &\cdots &f_m(x_n) \\ \end{bmatrix} </math> Med

<math>

\mathbf c =\begin{bmatrix} c_1 \\ c_2\\ \vdots \\ c_m\\ \end{bmatrix} ,\quad \mathbf y = \begin{bmatrix} y_1 \\ y_2\\ \vdots \\ y_n\\ \end{bmatrix} </math>

Fil:Offsets-datapoint-graf.svg

Avståndet ("felet") mellan datapunkter och approximerande kurva mäts i "vertikal" led och inte som punkternas vinkelräta avstånd till kurvan

kan ett linjärt ekvationssystem (vanligen överbestämt, normalt är n betydligt större än m) i m obekanta skrivas

<math>A\cdot \mathbf c = \mathbf y</math>

Att lösa detta ekvationssystem i minstakvadratmetodens mening är ekvivalent med att lösa normalekvationen

<math>A^T A\, \mathbf c = A^T\, \mathbf y</math>

där A^T är transponatet till A.

Om A och y har samma antal rader och om kolumnvektorerna i A är linjärt oberoende, har normalekvationen en entydig lösning c_min, för vilken gäller

<math>\| A\mathbf\, c - \mathbf y\|^2 \geq \| A\,\mathbf c_{min} - \mathbf y\|^2</math>

det vill säga, c_min är minimumpunkten till funktionen

<math>\mathbf c \rightarrow \|A\,\mathbf c - \mathbf y\|^2</math>

Det kvadratiska medelfelet beräknas som

<math>\epsilon = \| A\,\mathbf c_{min} - \mathbf y \|/\sqrt n</math>

Anpassning av polynom

För att anpassa ett polynom av grad m

<math> c_0 + c_1 x + \dots + c_m x^m</math>

till datamängden

<math>(x_1,\ y_1),\ (x_2,\ y_2),\,\dots,\ (x_n,\ y_n)</math>

sätts polynomets monom (med alla c_i = 1) med beräknade värden in som rader i designmatrisen

<math>A = \begin{bmatrix}1 &x_1^1 &\cdots &x_1^m\\

1 &x_2^1 &\cdots &x_2^m\\ \vdots &\vdots & \vdots & \vdots \\ 1 &x_n^1 &\cdots &x_n^m\end{bmatrix}</math> De sökta koefficienterna c och alla y-värden bildar kolumnvektorerna

<math>\mathbf{c}=\begin{bmatrix}c_0 \\ c_1 \\ \vdots \\ c_m \end{bmatrix},\quad

\mathbf{y}=\begin{bmatrix}y_1 \\ y_2 \\ \vdots \\ y_n \end{bmatrix}</math>

Därefter löses vanligen normalekvationen

<math>A^TA \cdot \mathbf{c} =A^T\cdot\mathbf{y}</math>

Fil:X22-multiple-grafs.png

Vilket gradtal skall väljas?

Val av polynomets grad

Givet värdet av datamängdens storlek, n, hur skall det approximerande polynomets grad m väljas? Grundantagandet är att m < n, eller åtminstone att datamängden med tillräcklig noggrannhet kan approximeras av ett sådant polynom. Om m ≥ n förbättras inte approximationen. Är m = n - 1 är lösningen exakt, men i detta fall förloras en vanligen önskvärd egenskap hos polynomet, nämligen förmågan att filtrera bort detaljer orsakade av mätfel och andra störningar (till exempel numeriska fel).

Normalekvationen

Fil:Normalekvationen.png

Ett vektorrum V spänns upp av A:s kolonnvektorer, i detta fall u och v. Vektorn b tillhör inte V varför ekvationssystemet A c = b saknar lösning. A c_min - b är ortogonal mot alla vektorer i V.
c_min kan ses som de koefficienter som minimerar "längden" av vektorn A c - b

Som en orientering beskrivs kortfattat en bakgrund till normalekvationen

<math>A^TA\,\mathbf c = A^T\mathbf b</math>

i form av ett specialfall (illustrerbart) med tre linjära ekvationer och två obekanta koefficienter. Antag att

<math>A = \begin{bmatrix}

u_1 & v_1\\ u_2 & v_2\\ u_3 &v_3\\ \end{bmatrix} \ =\ \begin{bmatrix} \mathbf u & \mathbf v\\ \end{bmatrix}, </math>

<math>

\mathbf c =\begin{bmatrix} c_1 \\ c_2\\ \end{bmatrix} ,\quad \mathbf b = \begin{bmatrix} b_1\\ b_2\\ b_3\\ \end{bmatrix}</math> och att kolonnvektorerna u och v i A spänner upp vektorrummet V (här, ett plan i R³). I allmänhet tillhör inte b vektorrummet V, varför ekvationen A c - b = 0 i allmänhet saknar lösning. Det är emellertid möjligt att söka en approximativ lösning, till exempel i minstakvadratmetodens mening, alltså en lösning till minimumproblemet

<math>\min_\mathbf c\,\|A\,\mathbf c - \mathbf b\|^2</math>

Detta minimum föreligger när A c - b är ortogonal mot vektorerna i V, det vill säga då skalärprodukterna av A c - b och varje vektor i V är noll. Men raderna i A:s transponat tillhör V och då matrisprodukten av A:s transponat och A c - b är definierad, ger detta

<math>A^T(A\,\mathbf c - \mathbf b) = 0 \quad\Rightarrow\quad A^TA\,\mathbf c = A^T\mathbf b</math>

och det sökta värdet på c, c_min, måste således satisfiera detta ekvationssystem. Om matrisen A^TA är inverterbar (om och endast om, kolonnerna i A är linjärt oberoende) är lösningen

<math>\mathbf c_{min} = (A^TA)^{-1}A^T\,\mathbf b</math>

och det går att visa att c_min uppfyller

<math>\| A\,\mathbf c - \mathbf b\|^2 \geq \| A\,\mathbf c_{min} - \mathbf b\|^2</math>

Dessa resultat är i huvudsak tillämpbara på allmänna rektangulära matriser A.

Lösningar om designmatrisens kolonner är ortogonala

Sök en lösning till ekvationen <math>A\cdot \mathbf c = \mathbf b</math> om

<math>A = \begin{bmatrix}1 &-12\\

1 &-4\\ 1 &2\\ 1 &14 \end{bmatrix} = \begin{bmatrix} \mathbf u & \mathbf v \end{bmatrix} ,\quad \mathbf b = \begin{bmatrix} -2\\ 4\\ 2\\ 12 \end{bmatrix}</math> Eftersom kolonnerna i A är ortogonala (<math>\mathbf u\cdot \mathbf v = 0</math>) ges den ortogonala projektionen av b på A:s kolonnrum av

<math>\hat{\mathbf b} =\frac{\mathbf b\cdot \mathbf u}{\mathbf u\cdot \mathbf u}\mathbf u + \frac{\mathbf b\cdot \mathbf v}{\mathbf v\cdot \mathbf v}\mathbf v =

4\mathbf u + \frac{1}{2}\mathbf v = \begin{bmatrix} 4\\ 4\\ 4\\ 4 \end{bmatrix} + \begin{bmatrix} -6\\ -2\\ 1\\ 7 \end{bmatrix} = \begin{bmatrix} -2 \\ 2 \\ 5 \\ 11 \end{bmatrix} \quad (1) </math> Då den ortogonala projektionen <math>{\mathbf \hat b}</math> är känd går det att lösa <math>A\cdot \mathbf{\hat c} = \mathbf{\hat b}</math>. Enligt (1) är <math>\mathbf{\hat c} = \begin{bmatrix} 4 \\ \cfrac{1}{2} \end{bmatrix} </math>, vilket i minstakvadratmetodens mening också är lösningen till <math>A\cdot \mathbf c = \mathbf b</math>.

Matriser där kolonnerna är ortogonala förekommer relativt ofta i problem inom linjär regression.

Exempel

Anpassning av en rät linje

Vilken rät linje

<math>y = c_1 x + c_0</math>

ger bästa anpassningen till mätserien

<math>(x_1,\ y_1),\ (x_2,\ y_2),\,...\,,\ (x_n,\ y_n)</math>

I detta fall blir designmatrisen

<math>A = \begin{bmatrix}1 &x_1\\ 1 &x_2\\ \vdots & \vdots \\1 &x_n\end{bmatrix}</math>

och y-värdena och de sökta koefficienterna placeras i

<math>\mathbf{y} = \begin{bmatrix}y_1 \\ y_2 \\ \vdots \\ y_n \end{bmatrix}, \quad

\mathbf{c} = \begin{bmatrix}c_0 \\ c_1 \end{bmatrix}</math> Därefter löses

<math>A^TA \cdot \mathbf{c} =A^T\cdot\mathbf{y}</math>

med avseende på c.

Anpassning av ett andragradspolynom

Fil:LeastSqrPoly.png

Exemplets polynomanpassning

Givet datapunkterna (1,10), (2,8), (3,11), (4,17), (5,24) söks de koefficienter till andragradspolynomet

<math>y = c_2x^2+c_1 x+c_0</math>

som enligt minstakvadratmetoden är bäst anpassade till observationerna.

Designmatrisen och vektorn för y-värdena är

<math>A = \begin{bmatrix}1 & 1 & 1 \\ 1 & 2 & 4 \\ 1 & 3 & 9 \\ 1 & 4 & 16 \\ 1 & 5 & 25 \end{bmatrix},\quad

\mathbf{y} = \begin{bmatrix}10 \\ 8 \\ 11 \\ 17 \\ 24 \end{bmatrix}</math>

<math>A^TA = \begin{bmatrix}5 &15 & 55 \\ 15 & 55 & 225 \\ 55 & 225 & 979 \end{bmatrix}, \quad A^T\mathbf{y}=\begin{bmatrix}70\\ 247 \\ 1013\end{bmatrix} </math>

Normalekvationen löses med avseende på c

<math>A^TA\begin{bmatrix}c_0 \\ c_1\\c_2\end{bmatrix} = A^T\mathbf{y} \quad \Rightarrow\quad\begin{matrix}c_0 = &13{,}4\\ c_1 = &-5{,}3 \\c_2 = &1{,}5 \end{matrix}</math>

och det anpassade andragradspolynomet är således

<math>y = 1{,}5x^2 -5{,}3x + 13{,}4 </math>

Jämförelse mellan observerade och minstakvadratanpassade y-värden.

x	uppmätt y	anpassat y	felet	felet i kvadrat
1	10	9,6	-0,4	0,16
2	8	8,8	0,8	0,64
3	11	11,0	0,0	0,00
4	17	16,2	-0,8	0,64
5	24	24,4	0,4	0,16
Summa:				1,60

Av alla möjliga andragradspolynom har inget en summa av felen i kvadrat som understiger 1,6.

Anpassning av ellips

Fil:LeastSqr-ellips.png

Den anpassade ellipsen

Kan datapunkterna (-9, 2), (-2, 5), (3, 6), (7, 4), (9, 1), (8, -4), (1, -5), (-4, -5), (-8, -3), (-9, -1) på ett meningsfullt sätt beskrivas av en ellips? Minstakvadratmetoden kan användas för att anpassa en ellips till datamängden. Ekvationen för en ellips är

<math>\frac{x^2}{a^2}+ \frac{y^2}{b^2}=1</math>

där a, b är ellipsaxlarnas längder.

De beräknade värdena för ellipsekvationens termer (med a och b = 1) sätts in i designmatrisens rader och värdena i ellipsekvationens högerled sätts in i kolumnvektorn b:

<math>A = \begin{bmatrix}x_1^2 & y_1^2 \\x_2^2 & y_2^2 \\ \vdots &\vdots \\ x_n^2 &y_n^2\end{bmatrix}\ = \begin{bmatrix}81 & 4 \\4 & 25 \\ \vdots &\vdots \\ 81 &1\end{bmatrix}, \quad \mathbf{b} = \begin{bmatrix}1 \\ 1 \\ \vdots \\ 1 \end{bmatrix}</math>

<math>A^TA = \begin{bmatrix}30630 & 3719 \\ 3719 & 3782\\ \end{bmatrix}, \quad A^T\mathbf{b}=\begin{bmatrix}450\\ 158\end{bmatrix} </math>

Normalekvationen löses

<math>A^TA\begin{bmatrix}c_1\\c_2\end{bmatrix}=A^T\mathbf{b}\quad\Rightarrow\quad\begin{matrix}c_1 = 0.010923 \\ c_2 = 0.031036\end{matrix}</math>

och därmed är

<math>a = \sqrt{\frac{1}{c_1}} = 9.5681</math>

<math>b = \sqrt{\frac{1}{c_2}} = 5.6764</math>

Anpassning av en yta

Fil:X22-3D-1.png

Den anpassade ytan

Anpassning av en yta i R³,

<math>z(x,y)=c_1\,x^3+c_2\,y^3+c_3\,x\,y</math>

till datapunkterna (x, y, z-koordinater i R³)

(2, 4, 33), (-1, 1, 2), (1, -3, 7), (4, 4, 88), (-2, -3, 26), (-3, 1, 13), (-1, -1, 4), (4, 1, 36)

Designmatrisen A konstrueras och datapunkternas z-värden placeras i kolonnvektorn z:

<math>A = \begin{bmatrix}x_1^3 & y_1^3 &x_1\,y_1 \\x_2^3 & y_2^3 &x_2\,y_2\\ \vdots &\vdots &\vdots \\ x_n^3 &y_n^3 &x_n\,y_n\end{bmatrix}\ = \begin{bmatrix} 8 & 64 &8 \\-1 & 1 &-1 \\ \vdots &\vdots &\vdots\\64 &1 &4\end{bmatrix}</math>

<math>\mathbf{z} = \begin{bmatrix}z_1 \\ z_2 \\ \vdots \\ z_n \end{bmatrix} = \begin{bmatrix}33 \\ 2 \\ \vdots \\ 36\end{bmatrix}</math>

Normalekvationen kan nu ställas upp och lösas:

<math>A^TA\begin{bmatrix}c_1\\c_2\\c_3 \end{bmatrix}=A^T\mathbf{z}\quad\Rightarrow\quad\begin{matrix}c_1 = &0.218245\\c_2 =&-0.033352\\c_3 = &4.241390\end{matrix}</math>