Linear least squares

Joshua R. Davis, Carleton College, Math 134, Fall 2024

Introduction

At its most basic level, linear algebra is about solving a linear system A x = b for the unknown vector x. But what if there is no solution? In a real-world situation, maybe we should still try to find something close to a solution. Maybe that “approximate solution” still has practical value.

Here’s the standard approach to such problems. Find the c that is closest to b, for which the system A x = c has a solution, and then find that solution x. This approach is called linear least squares. It is the heart of linear regression, which is an essential tool of statistics and machine learning.

Infrastructure

Run the huge cell of code below, but don’t worry about how it works. (Remember that, to run a cell, you click anywhere in the cell, and then press Shift-Enter.)

Clear[transformedCurvesPoints, transformer];
(* Helper function for transformer. *)
transformedCurvesPoints[bound_, curves_, points_, matrix_] :=
	Module[{objects},
		objects = {
			RGBColor[0.0, 0.0, 0.0],
			PointSize[Large],
			Point[Map[Function[point, Chop[point.Transpose[matrix]]], points]],
			AbsoluteThickness[2],
			Line[Map[Function[curve, Chop[curve.Transpose[matrix]]], curves]]};
		If[Length[matrix] == 2,
			Graphics[Join[{
				RGBColor[1.0, 0.0, 0.0],
				Line[{{bound, 0}, {0, 0}}],
				Line[{{-bound, 0}, {-bound/2, 0}}],
				RGBColor[0.0, 1.0, 0.0],
				Line[{{0, bound}, {0, 0}}],
				Line[{{0, -bound}, {0, -bound/2}}]},
				objects]],
			Graphics3D[Join[{
				RGBColor[1.0, 0.0, 0.0],
				Line[{{bound, 0, 0}, {0, 0, 0}}],
				Line[{{-bound, 0, 0}, {-bound/2, 0, 0}}],
				RGBColor[0.0, 1.0, 0.0],
				Line[{{0, bound, 0}, {0, 0, 0}}],
				Line[{{0, -bound, 0}, {0, -bound/2, 0}}],
				RGBColor[0.0, 0.0, 1.0],
				Line[{{0, 0, bound}, {0, 0, 0}}],
				Line[{{0, 0, -bound}, {0, 0, -bound/2}}]},
				objects]]]];
(* Returns a manipulator for a 2D or 3D linear transformation. *)
(* bound determines the length of the coordinate axes and often the view. *)
(* curves and points are some curves and points to be transformed. *)
(* matrix is a t-dependent matrix used to transform the curves and points. *)
transformer[bound_, curves_, points_, matrix_] :=
	Manipulate[transformedCurvesPoints[bound, curves, points, matrix], {t, 0, 1}];
(* Simple 2D and 3D geometry to transform. *)
box2D = {{{1, 1}, {-1, 1}, {-1, -1}, {1, -1}, {1, 1}}, 
{{1, 0}, {0, 0}}, {{0, 1}, {0, 0}}};
circlePoints[{a_, b_}, r_, n_] := 
Table[{r Cos[t] + a, r Sin[t] + b}, {t, 0, 2 Pi, 2 Pi / n}];
circle2D = circlePoints[{0, 0}, 1, 100];
rock2D = Append[box2D, circle2D];
box3D = {{{-1,-1,-1},{-1,-1,1}},{{-1,1,-1},{-1,1,1}},{{1,-1,-1},{1,-1,1}},{{1,1,-1},{1,1,1}},{{-1,-1,-1},{-1,1,-1}},{{-1,-1,1},{-1,1,1}},{{1,-1,-1},{1,1,-1}},{{1,-1,1},{1,1,1}},{{-1,-1,-1},{1,-1,-1}},{{-1,-1,1},{1,-1,1}},{{-1,1,-1},{1,1,-1}},{{-1,1,1},{1,1,1}},
{{0,0,0},{0,0,1}},{{0,1,0},{0,1,1}},{{1,0,0},{1,0,1}},{{0,0,0},{0,1,0}},{{0,0,1},{0,1,1}},{{1,0,0},{1,1,0}},{{0,0,0},{1,0,0}},{{0,0,1},{1,0,1}},{{0,1,0},{1,1,0}}};

Projection onto the column space of A

Enter any two “random” 3D vectors a1 and a2. (For now, just don’t make them scalar multiples of each other.) We’ll put them into the columns of a 2x3 matrix A.

a1 = {1, -1, 0};
a2 = {2, 0, 1};
a = Transpose[{a1, a2}];
MatrixForm[a]

If the 2x2 matrix A^T A is invertible --- as it usually is, in applications of this technique --- then we can compute the 3x3 matrix P = A (A^T A)^-1 A^T, and that matrix represents orthogonal projection onto the column space of A.

p = a.Inverse[Transpose[a].a].Transpose[a];
MatrixForm[p]

What is the following code checking?

p.a1
p.a2

Let’s visualize what the matrix P does to 3D space.

Clear[t];
MatrixForm[p]
transformer[3, box3D, {}, IdentityMatrix[3] + t (p - IdentityMatrix[3])]

How does this relate to the Introduction above? Well, when we can’t solve A x = b, we can (often) make P from A, and then solve A x = P b.

Go back and change a1 and a2, to see how they affect P. Moreover, try an example where a1 and a2 are scalar multiples of each other. See how the calculation breaks.

Fitting a polynomial: Math version

Suppose that we have a set of four data points (x1, y1), (x2, y2), (x3, y3), (x4, y4).

{x1, y1} = {2, 0};
{x2, y2} = {3, -1};
{x3, y3} = {-1, 0};
{x4, y4} = {5, 1};
ListPlot[{{x1, y1}, {x2, y2}, {x3, y3}, {x4, y4}}]

We would like to fit a curve to them. On a purely mathematical level, we should be able to fit a cubic polynomial, meaning a function of the form y = a x^3 + b x^2 + c x + d. We just have to find the as-yet-unknown coefficients a, b, c, d. Asking the polynomial to hit those four points amounts to four linear equations in the four unknowns.

Clear[a, b, c, d];
y1 == a x1^3 + b x1^2 + c x1 + d
y2 == a x2^3 + b x2^2 + c x2 + d
y3 == a x3^3 + b x3^2 + c x3 + d
y4 == a x4^3 + b x4^2 + c x4 + d

Let’s rewrite that system in matrix form. Ordinarily our notation would be something like A x = b. But in this case the matrix A is made from the xi, so we’re going to call it X. And the vector b is made from the yi, so we’re going to call it y. And the unknown vector x consists of a, b, c, d, so we’re going to call it a. So our matrix equation is X a = y, where X and y are known and a is unknown.

xx = Map[Function[x, {x^3, x^2, x, 1}], {x1, x2, x3, x4}];
MatrixForm[xx]
yy = {y1, y2, y3, y4};
MatrixForm[yy]

Is the matrix X invertible? Yes. So the solution is easy to compute.

sol = Inverse[xx].yy;
MatrixForm[sol]

Here’s a picture of how well we did.

Show[
ListPlot[{{x1, y1}, {x2, y2}, {x3, y3}, {x4, y4}}], 
Plot[sol[[1]] x^3 + sol[[2]] x^2 + sol[[3]] x + sol[[4]], 
{x, Min[{x1, x2, x3, x4}], Max[{x1, x2, x3, x4}]}]]

Try changing the data points and examining how they change the cubic fit. You might even try an example where all of the data points happen to lie along a line. You might even try an example where one data point is exactly above another. But don’t change how many data points there are, yet.

Fitting a polynomial: Stats version

Now suppose that we add a fifth data point (x5, y5) and repeat the key calculations. The calculations break, because X is 5x4 and hence definitely not invertible. Algebraically, it’s because there’s no reason to expect a solution to a system of five equations in four unknowns (except in special cases). Geometrically, it’s because there’s no reason to expect a cubic polynomial to pass through five specified points (except in special cases).

{x1, y1} = {2, 0};
{x2, y2} = {3, -1};
{x3, y3} = {-1, 0};
{x4, y4} = {5, 1};
{x5, y5} = {4, 1};
xx = Map[Function[x, {x^3, x^2, x, 1}], {x1, x2, x3, x4, x5}];
MatrixForm[xx]
yy = {y1, y2, y3, y4, y5};
MatrixForm[yy]
sol = Inverse[xx].yy;
MatrixForm[sol]

At this point, a mathematician might give up, but a statistician would not. How can the statistician proceed? Here’s the key: The 4x4 matrix X^T X is invertible.

MatrixForm[Transpose[xx].xx]
MatrixForm[Inverse[Transpose[xx].xx]]

Therefore P = X (X^T X)^-1 X^T can be computed, and it represents the projection onto the column space of X. The linear least squares “solution” to X a = y is the solution to X a = P y. This solution is a = (X^T X)^-1 X^T y.

p = xx.Inverse[Transpose[xx].xx].Transpose[xx];
MatrixForm[p]
sol = Inverse[Transpose[xx].xx].Transpose[xx].yy;
MatrixForm[sol]

Examine how well we did in fitting the five data points.

Show[
ListPlot[{{x1, y1}, {x2, y2}, {x3, y3}, {x4, y4}, {x5, y5}}], 
Plot[sol[[1]] x^3 + sol[[2]] x^2 + sol[[3]] x + sol[[4]], 
{x, Min[{x1, x2, x3, x4, x5}], Max[{x1, x2, x3, x4, x5}]}]]

Here’s an exercise for you: Can you rig the five data points so that X^T X is not invertible and therefore the least squares calculation breaks? If this task seems daunting, then return to earlier sections of this notebook, to see how earlier calculations were broken.

More data

Enter any cubic polynomial you like. (The pictures work better if you use small coefficients.) I will make a synthetic data set of 100 points that “noisily” follow that polynomial.

Clear[x];
polyn = x^3 - 4 x^2 + 2 x + 1;
xs = RandomVariate[UniformDistribution[{-3, 3}], 100];
ys = RandomVariate[NormalDistribution[0, 5], 100];
ys = ys + Map[Function[pred, (polyn /. x -> pred)], xs];
data = Transpose[{xs, ys}];
ListPlot[data]

Now our X and y matrices are big: 100x4 and 100x1, respectively.

xx = Map[Function[x, {x^3, x^2, x, 1}], xs];
MatrixForm[xx]
yy = ys;
MatrixForm[yy]

That X matrix is nowhere close to invertible. But X^T X is still 4x4 and invertible.

MatrixForm[Transpose[xx].xx]
MatrixForm[Inverse[Transpose[xx].xx]]

The projection matrix P is now 100x100. The solution vector a is still 4x1.

p = xx.Inverse[Transpose[xx].xx].Transpose[xx];
MatrixForm[p]
sol = Inverse[Transpose[xx].xx].Transpose[xx].yy;
MatrixForm[sol]

Let’s see how well we did in fitting the 100 data points.

Show[
ListPlot[data], 
Plot[sol[[1]] x^3 + sol[[2]] x^2 + sol[[3]] x + sol[[4]], {x, Min[xs], Max[xs]}]]

The word “linear”

Thus far we’ve been fitting cubic polynomials all the time. But linear least squares can be used to fit many kinds of curves (and even higher-dimensional models, such as planes and surfaces). People really like fitting lines to data, because lines are easy to interpret. Here’s a line fit to the data set from the preceding section.

xx = Map[Function[x, {x, 1}], xs];
yy = ys;
sol = Inverse[Transpose[xx].xx].Transpose[xx].yy;
MatrixForm[sol]
Show[
ListPlot[data], 
Plot[sol[[1]] x + sol[[2]], {x, Min[xs], Max[xs]}]]

A statistician would tell you that the line is not a good fit, because the data set has a visible curving tendency that the line doesn’t capture. They might suggest that you try a quadratic curve y = a x^2 + b x + c or a cubic curve or a sinusoidal curve or an exponential curve or some other kind of curve. How do we know, which kind of curve is the right one for the data, if we don’t know how the data were constructed? That’s tough. Learning how to use linear least squares well takes a lot of practice. Anyway, that’s not my point.

My point is this: Sometimes students come away with the misconception that linear least squares is called “linear” because it’s about fitting lines to data. No, that is not right. In most of this notebook, we’ve been using linear least squares to fit non-linear functions to data, right? Linear least squares does not require that the y-values depend linearly on the x-values. Rather, it requires that the y-values depend linearly on the coefficients a.

An exercise

This data set is explained in Discussion Problem 13. The q-values are quantities of an element observed at times t.

ts = {0, 1, 2, 3, 4, 5, 6};
qs = {1.284,1.160,1.119,0.860, 0.858,0.787,0.692};
data = Transpose[{ts, qs}];
ListPlot[data, AxesOrigin->{0, 0}]

Although the graph looks fairly linear, theory predicts that it should look like exponential decay. That is, there should be numbers a and b such that q = a e^(b t). If we set z = log(q) and c = log(a), where log denotes the natural or base-e logarithm, then the relationship becomes z = c + b t. Here are the z-values and the q-values re-computed from the z-values as a check.

zs = Log[qs]
E^zs

Your job is: By copying, pasting, and adjusting the code from the previous section of this notebook, use linear regression to discover c and b. Then compute a = e^c. Once you know a and b, you can type them into the code cell below, to see how good your fit is. (The values a = 1 and b = 0 are placeholders; they are not correct.)

a = 1.0;
b = 0.0;
Show[
ListPlot[data], 
Plot[a E^(b t), {t, Min[ts], Max[ts]}]]

In practice

In practice, if you want to fit models to data using linear least squares, don’t go through the calculations manually, as we did above. Why? The obvious reason is that this fitting problem is so common that there must be a method built into Mathematica. The less obvious, but more important, reason is that numerical calculations on a computer are imprecise and error-prone, and there are ways to organize the calculations that are more error-resistant than what we did above. Maybe we’ll talk about them later in this course.

Anyway, here is how you might actually fit the cubic polynomial to the data above, if your job depended on it. But notice that we get almost exactly the same answer.

lm = LinearModelFit[data, {1, x, x^2, x^3}, x]
Show[
ListPlot[data], 
Plot[lm[x], {x, Min[xs], Max[xs]}]]