Strang's Linear Algebra: least squares

Strang has an excellent section on least squares. I'm going to work through an example here as a way of building my understanding (homework). We take the simplest possible case of three time points equally spaced at t = 1, 2, 3. At these times we observe data y = 2, 2, 4 (to make it slightly different than his example). The three points are not co-linear, as you might guess and can see from the plot. The numpy function polyfit will fit a line to this data (that's what I plotted above), but we're going to find it ourselves.

The line that "best fits" the data has the equation y = C + Dt. We need to find C and D. We use calculus first, then linear algebra. The magnitude of the "error" at times 1, 2 and 3 is the difference between the value on the line and the value we actually observed.

e1 = C + 1D - 2 e2 = C + 2D - 2 e3 = C + 3D - 4

We calculate the sum of the squares of the errors, which we wish to minimize. One reason for the square is that the minimization involves the derivative, and the derivative of the square function is a linear function.

SumSq = (C+D-2)2 + (C+2D-2)2 + (C+3D-4)2

This is a function of two variables (C and D), so we need to take both partial derivatives and set them both equal to zero, (yielding two equations in two unknowns).

∂SumSq/∂C = 2(C+D-2) + 2(C+2D-2) + 2(C+3D-4) = 0 ∂SumSq/∂D = 2(C+D-2) + 2*2(C+2D-2) + 3*2(C+3D-4) = 0

The factor of 2 in each term goes away. Notice that the first equation says that the sum of the errors is equal to zero.

Solve for C in the first eqn:

C + D - 2 + C + 2D - 2 + C + 3D - 4 3C + 6D = 8 3C = 8 - 6D

Substitute C into the second equation:

C + D - 2 + 2C + 4D - 4 + 3C + 9D - 12 = 0 6C + 14D = 18 14D = 18 - 6C 14D = 18 - 16 + 12D 2D = 2 D = 1 C = 2/3

In Bolstad (p. 237) the equation for the slope is given as:

B = mean(xy) - mean(x)*mean(y) / mean(x2) - mean(x)2 mean(xy) = mean(1*2 + 2*2 + 3*4) = 18/3 = 6 mean(x)*mean(y) = 2 * mean(2,2,4) = 16/3 mean(x2) = mean(1,4,9) = 14/3 mean(x)2 = 4 B = 18/3 - 16/3 = 1 ----------- 14/3 - 12/3

and the intercept is

Ao = mean(y) - B*mean(x) = 8/3 - 2 = 2/3

Looks good. Let's do some linear algebra! We go back to the equation for the line: C + Dt. We have one equation for each observation, which is written like so:

[ 1 1 ] [ C ] = [ 2 ] [ 1 2 ] [ D ] [ 2 ] [ 1 3 ] [ 4 ]

Ax = b has no solution. We seek the projection of b onto A.

Using the formula from last time, (and Python to do the calculation, below) we obtain:

P = [ 5/6 1/3 -1/6 ] [ 1/3 1/3 1/3 ] [ -1/6 1/3 5/6 ] p = P b = [ 5/3 ] [ 8/3 ] [ 11/3 ]

We can check that the error e = b - p

[ 2 ] - [ 5/3 ] = [ 1/3 ] [ 2 ] - [ 8/3 ] = [ -2/3 ] [ 4 ] - [ 11/3 ] = [ 1/3 ]

is perpendicular to p:

e•p = 5/9 - 16/9 + 11/9 = 0

import numpy as np import matplotlib.pyplot as plt t = np.array([1,2,3]) b = np.array([2,2,4]) plt.scatter(t,b,s=300,color='r') ar,by = np.polyfit(t,b,1) def f(x): return by + ar*x T = np.arange(0,3.2,0.1) plt.plot(T,f(T),lw=4,color='k') plt.text(0.2,3.5,fontsize=24, s='y-intercept =' + str(round(by,3))) plt.text(0.2,3.0,fontsize=24, s='slope =' + str(round(ar,3))) ax = plt.axes() ax.set_xlim(0,3.2) ax.set_ylim(0,4.2) plt.savefig('example.png') col1 = np.array([1,1,1]) t.shape = b.shape = col1.shape = (3,1) A = np.hstack([col1,t]) print A def project(A): AT = A.transpose() temp = np.linalg.inv(np.dot(AT,A)) return np.dot(A,np.dot(temp,AT)) P = project(A) print P print b print np.dot(P,P) print np.dot(P,b)