Least Squares

We'll start with some data, which we want to fit to a model. Note how vectors and matrices are defined in Mathematica. A vector is a list enclosed by curly braces. And a matrix is a list of lists.

data = {{0,0},{Pi/3, Sqrt[3]/2},{Pi/2, 1},{5 Pi/6,1/2},{Pi,0}}

datapoints = ListPlot[data,PlotStyle->PointSize[0.02]]

Will a quadratic function fit these data? We'll assume such a model: y = a + b x + c (x^2).

A = {{1,0,0},{1, Pi/3, (Pi/3)^2},{1, Pi/2, (Pi/2)^2},{1, 5Pi/6, (5Pi/6)^2},{1, Pi, Pi^2}}; A // MatrixForm

A.{a,b,c} // MatrixForm

bvec = {0, Sqrt[3]/2, 1, 1/2, 0}; bvec//MatrixForm

If we naively try to solve the linear system Ax = b, we get

LinearSolve[A,bvec]

The Least Squares solution is

x = Inverse[Transpose[A].A].Transpose[A].bvec

We prefer a numerical solution.

xx = N[x]

The model is thus

model = xx . {1, z, z^2}

Plotting the data with the model gives

plot = Plot[model,{z,0,Pi}]

Show[datapoints,plot]

Of course, we know a better fit for these data.

Show[datapoints,Plot[{Sin[z],model},{z,0,Pi}]]

What about trying a cubic function instead of a quadratic? y = a + b x + c (x^2) + d (x^3). The A matrix now becomes

A = {{1,0,0,0},{1, Pi/3, (Pi/3)^2, (Pi/3)^3},{1, Pi/2, (Pi/2)^2, (Pi/2)^3},{1, 5Pi/6, (5Pi/6)^2, (5 Pi/6)^3},{1, Pi, Pi^2, Pi^3}}; A // MatrixForm

That's the only thing we need to change. Now just plug in the least-squares formula:

x = Inverse[Transpose[A].A].Transpose[A].bvec

N[x]

newmodel = N[x] . {1, z, z^2, z^3}

plot = Plot[newmodel,{z,-30,30}]

plot = Plot[newmodel,{z,0,Pi}]

Show[datapoints,plot]

Of course, we know a better fit for these data.

Show[datapoints,Plot[{Sin[z],newmodel},{z,0,Pi}]]

Consider the height and weight data collected in class last week.

height = {68.,76,72,68,70,74,62,66,67,72,66,68,68,70,72}

weight= {152.,217,170,125,146,175,115,116,140,183,155,133,148,149,160}

Transpose[{height,weight}]

scatter =ListPlot[Transpose[{height,weight}], PlotStyle->PointSize[0.02]]

To obtain the least squares line, we consider "solving" the linear equations of the form WEIGHT = a + b (HEIGHT) . That is, we have a linear system of the form Ax = P, where x = [ a, b ] is the vector of unknown parameters for the line we're trying to find. P is the WEIGHT vector. And the coefficient matrix A is

A = Transpose[{height^0.,height}]; A//MatrixForm

A.{a,b} //  MatrixForm

Using the least squares formula gives

x = Inverse[Transpose[A].A].Transpose[A].weight

equation = x.{1,z}

The model is telling us that for every inch increase in height, weight goes up by about 6.5 pounds.

leastsquaresline = Plot[equation,{z,60,80}]

Show[scatter,leastsquaresline]

Perhaps a cubic polynomial would give a better fit? Consider the cubic WEIGHT = a + b (HEIGHT) + c (HEIGHT)^2 + d ( HEIGHT)^3.

Our new coefficient matrix is

A = Transpose[{height^0,height, height^2,height^3}]

Plug and chug away . . .

Inverse[Transpose[A].A].Transpose[A].weight

Oh, oh numerical problems. . . . We'll try to get around it by considering the "normalized" height and pounds vectors obtained by subtracting off the mean. Thus our two vectors will be "centered" around 0 and we shouldn't expect exceptionally large or small values from squaring and cubing.

newheight = height-Mean[height]; newweight = weight-Mean[weight];

newscatter =ListPlot[Transpose[{newheight,newweight}], PlotStyle->PointSize[0.02],PlotRange->All]

A = Transpose[{newheight^0,newheight, newheight^2,newheight^3}]

Plug and chug away . . .

x =Inverse[Transpose[A].A].Transpose[A].newweight

Clear[z]

equation = x.{1,z,z^2, z^3}

leastsquarescubic = Plot[equation,{z,-8,8},PlotRange->All]

Show[newscatter,leastsquarescubic]

A = Transpose[{ newheight^0 , newheight,Cos[newheight], Sin[newheight], Cos[2 newheight], Sin[2 newheight]}]

x = Inverse[Transpose[A].A].Transpose[A].newweight

Clear[z]

equation = x.{1,z,Cos[z],Sin[z],Cos[2z],Sin[2z]}

leastsquares = Plot[equation,{z,-8,8},PlotRange->All]

Show[newscatter,leastsquares]

Find the Orthogonal Projection onto the plane x + 2y - z = 0. Check that the plane is the span of the vectors [-2, 1, 0] and [1, 0, 1].

A = {{-2,1},{1,0},{0,1}}; A//MatrixForm

OG = A.Inverse[Transpose[A].A].Transpose[A]; OG //MatrixForm

Recall that Projection matrixes have the property that they are equal to their own squares.

OG.OG // MatrixForm

The above two vectors v1 and v2 are a basis for the plane, but they are not an orthonormal basis. The Mathematica command Orthogonalize will return an orthonormal basis from a given set of vectors. (It uses the Gram-Schmidt procedure which we skipped in Section 5.2)

{v1, v2} =Orthogonalize[{{-2,1,0},{1,0,1}}]

Check that these vectors are orthonormal. The "norm" of a vector is the length of the vector.

Norm[v1]

Norm[v2]

v1.v2

Now check that the two vectors are actually on the plane.

v1.{1,2,-1}

v2.{1,2,-1}

With an orthonormal basis we get a simpler formula for the matrix of the orthogonal projection. Fact 5.3.10 gives matrix as (QQ^T), where Q is the matrix whose columns consist of the orthonormal basis vectors.

Q = Transpose[{v1,v2}]; Q //MatrixForm

Q.Transpose[Q] // MatrixForm