Least Squares
Math 272  Prof. Leise

Email completed labs to tleise@amherst.edu by 4pm Monday.

Today we will use least squares to "solve" inconsistent systems by finding the vector x that comes closest to solving Ax=b.  When fitting a line to data, least squares is often called linear regression (the basic idea is to minimize the sum of squares error between the data and the line).

To fit a line to a set of datapoints, we plug each point into y=ax+c to generate a set of linear equations in a and c.  In matrix form Ax=b, the matrix A will have the x's forming the first column and 1's forming the second column, while the vector b will consist of the y's.  The vector x={a,c} contains the slope and intercept of the line.  Unless the data are perfectly linear, the resulting system Ax=b will be inconsistent and we instead seek an approximate solution that minimizes the sum of squares error.  In linear algebra lingo, we want to solve the normal equation to find the least squares solution.

Recall that the normal equation for Ax=b is (A^T)Ax=(A^T)b.  You may use the LinearSolve command to solve the normal equation.

A={{1,1},{1,-1}};
b={1,-3};
LinearSolve[Transpose[A].A,Transpose[A].b]

Exercise 1.  Use the normal equation to an appropriate linear system Ax=b to find the best-fit line to the following set of data:

Gestation period in days  	Longevity in years 	(data from 1993 World Almanac and Book of Facts)
Black bear	219		18
Grizzly bear	225		25
Polar bear	240		20
Leopard	  98		12
Lion		100		15
Puma		  90		12
Tiger		105		16

Plot the data and the best-fit line using the command below (you need to enter the correct coefficients for the line marked in red).

data={{219,18},{225,25},{240,20},{98,12},{100,15},{90,12},{105,16}};

Show[Plot[0.5x+1,{x,0,250}],ListPlot[data,PlotStyle->PointSize[0.02]]]

To avoid lots of tedious data entry, we can use the Table command to generate a new list of points, a matrix, or a vector.  For example, the commands below generate the matrix A whose first column are the x's of the datapoints and whose second column consists of 1's, and the vector b consisting of the y's of the datapoints.  Note that data[[k,1]] gives the x-component of the kth datapoint and data[[k,2]] gives the y-component of the kth datapoint.

NumPoints=7;
A=Table[{data[[k,1]],1},{k,1,NumPoints}] 
b=Table[data[[k,2]],{k,1,NumPoints}]

Exercise 2.  Plot the following set of planetary data.  The relationship between distance and period of revolution is not linear, but it is a power.   Once you have the data entered as a list of points, use the Table trick to generate the new lists of data with powers of the distance: powerdata=Table[{data[[k,1]]^p,data[[k,2]]},{k,1,9}].  Then experiment with taking powers of the distance, for example, p=0.5, 1, 1.5, 2,..., until you find one that makes the ListPlot appear linear.

Average distance from Sun  	Period of revolution around Sun 
(in millions of miles)		(in Earth days)
Mercury	    36			88
Venus		    67			225
Earth		    93			365
Mars		  142			687
Jupiter		  484			4332
Saturn		  887			10760
Uranus		1765			30684
Neptune	2791			60188
Pluto		3654			90467

Exercise 3.  Use the normal equation for an appropriate system Ax=b to find the best-fit relation period=a(distance^p)+c, using the power p found in Exercise 2.   Plot this relation on the same graph as the datapoints.

The last exercise is from Visual Linear Algebra by Eugene Herman, Mike Pepe, and Eric Schulz.

Exercise 4. Suppose the proposed path of a natural gas pipeline runs near a fault plane that has caused some small earthquakes. The figure below shows a map of the area (a square, 500 km. on each side) and the proposed path of the pipeline.  The yellow dot marks the center of the blue square with coordinates (0,0,0).  In this coordinate system, the proposed pipeline passes through points with coordinates (0,-75,0) and (250,-125,0).

({{}})

Even a small earthquake could damage the pipeline, if the fault plane ruptures the surface and if it intersects the ground (z=0) on a line (called the "fault trace") that cuts through the pipeline. To determine the approximate location of the fault plane, geologists have recorded the locations of 32 seismic events (i.e., tiny movements picked up by sensitive instruments) along the subterranean (z<0) fault plane. The data below are the coordinates (x,y,z) of these seismic events.  We want to fit a plane z=ax+by+c to this data.  This plane will approximate the fault and we can then find the intersection of this fault plane with the surface z=0 to determine whether the fault plane is likely to intersect the pipeline (a very worrisome scenario).

data={{39.757,172.063,-57.8}, {15.208,8.778,-40.8}, {76.684,74.344,-45.5}, {36.412,16.465,-43.2}, {10.187,47.710,-35.0}, {47.028,66.459,-46.6}, {55.266,69.793,-47.9}, {65.724,79.812,-46.0}, {54.524,58.679,-47.9}, {73.135,96.507,-56.2}, {-20.092,30.219,-30.8}, {45.631,163.161,-57.3}, {75.724,10.939,-46.0}, {28.84,199.888,-59.6}, {40.205,24.236,-44.9}, {43.280,64.241,-48.5}, {-34.466,-25.164,-35.1}, {60.461,99.808,-50.2}, {65.669,102.041,-50.3}, {54.530,46.453,-46.0}, {73.218,74.278,-54.1}, {-70.042,93.306,-36.1}, {0,0,-35.20}, {-69.320,19.910,-42.8}, {77.261,5.389,-49.5}, {-83.101,109.166,-36.0}, {30.240,-19.078,-39.6}, {-9.504,33.434,-34.4}, {105.594,52.260,-48.1}, {50.785,153.154,-54.1}, {45.545,80.909,-43.3}, {5.243,-7.822,-38.7}};

(a) Use our least-squares method to find the plane that best fits the data.  (Use the Table trick again to generate A and b.)
(b) Use the code below to graph the proposed pipeline and the fault trace.  You will need to enter the correct expression for the fault plane (the red coefficients).
(c) Should there be concern about the pipeline based on this (oversimplifed) analysis? (We are, of course, making some pretty big assumptions, e.g., assuming that the fault plane is essentially linear.)

Clear[x,y,z]
fault = Plot3D[ -.1x-.1y-20,{x,-250,250},{y,-250,250},ColorFunction->Function[z,Blue],Mesh->False];
ground = Plot3D[ 0,{x,-250,250},{y,-250,250},ColorFunction->Function[z,Yellow],Mesh->False];
pts = Show[Graphics3D[{Red, AbsolutePointSize[6], Map[Point[#]&,data]}]];
pipe = 
ParametricPlot3D[{x,-.2x-75,0},{x,-250,250},PlotStyle->{Thickness[.02],Red}];
Show[{fault,ground,pts,pipe},AspectRatio->1]