Least Squares

Email completed labs to npflueger@amherst.edu by Wednesday night. Remember to CC all group members on the email. Also, include the last names of all group members in the name of the file, e.g. CrickWatson-lab2.nb.

Syntax interlude: the LinearSolve command

We’ve previously seen how to use the Solve and RowReduce commands to solve linear systems of equations. There is another command that can be used for this task, which is often convenenient when you are given the linear system as a matrix equation.

Suppose that you know the matrix A and vector b. Hit shift-enter on the line below to set these variables. Note: the semicolon at the end of each of these lines instructs Mathematica to execute the command without printing any output to the screen.

A = {{1,1,1},{1,2,4},{1,3,9}};
b = {2,5,10};

To solve the system Ax = b, the LinearSolve command can be used as follows. Hit shift-enter to see the result.

LinearSolve[A,b]

The normal equation and least-squares solutions

Today we will use a method called least squares to "solve" inconsistent systems by finding the vector x that comes closest to solving Ax=b. This is useful in situations where you know that a system will not have a solution, but you want the next-best thing.

Terminology: the normal equation for a (possibly inconsistent) linear system Ax=b is (A^T)Ax=(A^T)b.  The purpose of the normal equation is to provide a consistent linear system, whose solution are “as close as possible” to solutions of Ax = b (in a sense that we will make precise in class in a couple weeks). In other words, you may assume the following.

Fact: For any choice of A and b, the normal equation is consistent. Its solutions are “as close as possible” to solutions of the system A x = b, in a specific precise sense. (Take this fact for granted for now).

A common situation in which we might want an approximate solution to an inconsistent system is when we wish to fit a line to some data in the plane. When least-squares is applied to this problem, it is often called linear regression.

For example, suppose that we wish to fit a line to the following three points: (0,1), (1,5), and (2,11). To fit a line means that we want an equation y = u x + v whose graph passes through all three points. Viewing u and v as variables, we are, in effect, trying to solve this linear system of equations, for the variables u and v:
	0 u + v = 1
	1 u + v = 5
	2 u + v = 11
As a matrix equation, we are trying to solve A x = b, where A and b are as follows. Hit shift-enter below to see that this system is not consistent.

A={{0,1},{1,1},{2,1}};
b={1,5,11};
LinearSolve[A,b]

We can instead solve the normal equation, to get the least-squares solution, using the syntax below.

LinearSolve[ Transpose[A] . A, Transpose[A] . b]

Syntax comment: To find the transpose of a matrix A, write Transpose[A]. To multiply two matrices, use a period, as shown above.

Let’s consider an example with more data points. The table below shows the relationship between gestation period and longevity, for various mammal species. (Apologies if the indentation appears off on your screen; this issue seems to arise when a file is opened in a different version of Mathematica than it was created in.)

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

These data can be summarized in the following matrix.

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

To plot these data, use the following syntax. (Hit shift-enter in the line above to load “data” as a variable, then in the line below to make the plot using “data.”)

Show[ListPlot[data,PlotStyle->PointSize[0.02]]]

Exercise 1. Determine a matrix A and a vector b so that a solution to A x = b would amount to finding the coefficients of a line y = u x + v passing through these data points (follow the same logic as in the previous section). Note that, as in the example above, the second column of A should be all 1’s. Verify, using a method of your choice, that the equation A x = b has no solutions. See the syntax interlude after exercise 2 if you want to see a way to generate A and b a little faster than by simply typing them in.

Exercise 2.  Use the normal equation to find an approximate solution to A x = b. 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).

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

Syntax interlude: the Table command

In exercise 1 and 2, it isn’t too hard to type in the matrix A and the vector b by hand. In general, however, it could become tedious. Fortunately, there’s some handy Mathematica syntax that avoids this work: the Table command.

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. The way this syntax works is: the first argument to Table tells Mathematica what the kth row should look like. The second argument tells Mathematica which values of k to use to generate these rows (in this case, the numbers from 1 to NumPoints).

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

Exercise 3.  Plot the following set of planetary data. Use the ListPlot syntax from above.  For your convenience, you can write these numbers to the matrix “data” by hitting shift-enter in the line below the table. 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

data = {{36,88},{67,225},{93,365},{142,687},{484,4332},{887,10760},{1765,30684},{2791,60188},{3654,90467}};

Exercise 4.  Use the normal equation for an appropriate system Ax=b to find the best-fit relation period=u(distance^p)+v, using the power p found in Exercise 2.   Plot this relation on the same graph as the datapoints. You may see a warning about a badly conditioned matrix; you can ignore this warning.

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

Exercise 5. 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=ux+vy+w to this data, where u,v,w are the coefficients that we are hoping to solve for.  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.) Note: because we are now solving for three variables (u,v,w) instead of two, the matrix A should now have 3 columns instead of 2.

(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). If you have found the coefficents u,v,w correctly, the blue plan in this figure should very nearly pass through all of the red dots.

Clear[x,y,z]
fault = Plot3D[ -.1x-.1y-10,{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]

(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.)