A Second Course in Linear Algebra

Subsection 4.1.1 Theory of Least Squares

Solving a linear system of equations is a fundamental use of the tools of linear algebra. You know from introductory linear algebra that a linear system may have no solution. In an applied situation there could be many reasons for this, and it begs the question: what to do next?

We often construct mathematical models of practical situations, frequently in an effort to measure various parameters of the model. Suppose we think that interest rates \(R\text{,}\) as measured by the rate on one-year government bonds, are a linear function of construction activity, \(C\text{,}\) as measured by the number of permits issued for new construction in the last 30 days. So the “hotter” the construction market, the greater the demand for loans, so the cost of money (the interest rate) is greater. With a good model, we might be able to predict interest rates by examining public records for changes in the number of construction permits issued.

So we have a mathematical model

where we do not know \(a\) or \(b\text{,}\) the parameters of the model. But we would like to know the values of these parameters. Or at least have good estimates, since we understand that our model is an extremely simple representation of a much more complicated situation. So we collect data by obtaining monthly records of the interest rate and construction permits for the past sixty years, where maybe we report permits on a per capita basis to account for changes in population over many decades. Now we have \(720\) pairs of permits issued and interest rates. We can substitute each pair into our model and we get a linear equation in the parameters \(a\) and \(b\text{.}\) Two such equations would be likely to have a unique solution, however if we consider all \(720\) equations there is unlikely to be a solution. Why do we say that? Imagine the linear system we would obtain from all these equations. The \(720\times 2\) coefficient matrix will have the values of \(C\) in the first column, and the second column will be all \(1\)'s. The vector of constants will be the values of \(R\text{.}\) If you form the augmented matrix of this system and row-reduce, you will be virtually guaranteed to get a pivot column in the third column and thus have no solution (Theorem RCLS).

The lack of a solution for our model can be expressed more clearly by saying that vector of constants (the values of \(R\)) is not in the column space of the coefficient matrix. Or, the vector of \(R\) values is not a linear combination of the vector of \(C\) values and the vector of all ones. If it were, then the scalars in the linear combination would be our solution for \(a\) and \(b\text{.}\) We will temporarily leave our economic forecasting example to be more general, but we will return to the example soon.

Consider the general problem of solving \(A\vect{x}=\vect{b}\) when there is no solution. We desire an \(\vect{x}\) so that \(A\vect{x}\) equals \(\vect{b}\text{,}\) but we will settle for an \(\vect{x}\) so that \(A\vect{x}\) is very close to \(\vect{b}\text{.}\) What do we mean by close? One measure would be that we want the difference \(A\vect{x}-\vect{b}\) to be a short vector, as short as possible. In other words we want to minimize the norm, \(\norm{A\vect{x}-\vect{b}}\text{.}\) We will return to this minimization idea shortly, but let's ask informally what direction this short vector might have?

We have a subspace, the column space of \(A\text{,}\) represented by all possibilities for \(A\vect{x}\) as we course over all values of \(\vect{x}\text{.}\) We have a vector \(\vect{b}\) that is not in the subspace. Which element in the subspace is closest to \(\vect{b}\text{?}\) Consider the situation in two dimensions, where the subspace is a line through the origin, and the point lies off the line. From the point off the line, which of the infinitely many points on the line is closest, and what direction would you take to get there? In three dimensions, consider a subspace that is a plane containing the origin, and the point lies off the plane. From the point off the plane, which of the infinitely many points on the plane is closest, and what direction would you take to get there? You should have answered “in a perpendicular direction” and “in the direction of the normal vector to the plane”. Move orthogonally to the subspace. More exactly, in the direction of an element of the orthogonal complement of the subspace.

More carefully, let \(\hat{\vect{b}}\) be the “closest” vector to \(\vect{b}\) that is in the column space of \(A\text{,}\) and define \(\hat{\vect{x}}\) by the scalars in the linear combination of the columns of \(A\) that yields \(\hat{\vect{b}}\text{.}\) That is, \(A\hat{\vect{x}} = \hat{\vect{b}}\text{.}\) We desire

so

This linear system is called the normal equations, due to the role of orthogonality in its derivation. Several good things have happened. Primarily, \(\adjoint{A}A\) is a square matrix, is positive semi-definite, and if \(A\) has full column rank, then \(\adjoint{A}A\) is nonsingular. With a nonsingular coefficient matrix, we have the unique solution

The difference between \(\vect{b}\) and \(\hat{\vect{b}}\text{,}\) \(\vect{r}=\vect{b}-\hat{\vect{b}}\text{,}\) is known as the residual, and its length (the distance in the orthogonal direction) can be used as a measure of the approximation that is inherent in taking \(\hat{\vect{x}}\) as a solution to the system.

We consider two very simple examples of using the normal equations, in cases where our geometric intuition is useful.

Let's do it again, staying simple, but not trivially so.

What does a solution to the normal equations look like for our economic forecasting model? As discussed above \(A\) is a \(720\times 2\) matrix, where the first column has the numbers of construction permits, \(c_i\) for \(1\leq i\leq 720\text{,}\) and the second column is all ones. The vector \(\vect{b}\) contains the \(720\) values of the interest rate, \(r_i\text{,}\) \(1\leq i\leq 720\text{.}\) So

Then,

The expressions above can be cleaned up some, but the point is to see that we have an expression for each of \(a\) and \(b\text{,}\) that depends soley on the \(720\) pairs of data points \((c_i, r_i)\text{,}\) \(1\leq i\leq 720\text{.}\) These expressions may look familiar as the most basic case of linear regression. Exercise Checkpoint 4.1.3 asks you to derive these expressions, and in the more typical order. With estimates in hand, you can now consult the number of construction permits issued and form a prediction of interest rates.

So far we have not stated any theorems, to say nothing of proving anything. Moving in an orthogonal direction feels like a good idea, but is it really best? Here is a theorem that suggests it is best according to one natural criteria.

There are other measures of distance than the norm (we are using what would be called the \(2\)-norm, in order to differentiate it from other norms). But “Euclidean distance”, as we have used here, is often the most tractable and the results are often the most useful. In statistics, it is assumed that deviations in data, deviations from what the model would predict, are due to quantities that vary according to probability distributions. With this assumption, estimators, such as our least-squares estimates, become complicated functions of these random quantities. With the \(2\)-norm, and reasonable assumptions about the probability distributions for the data, the resulting least-squares estimators have probability distributions that are well-understood, thus making it possible to understand their properties and behavior.

Subsection 4.1.2 Goodness of Fit

We can associate a single number as a measure of how “good” our solution \(\hat{\vect{b}}\) really is. First, realize that the simplest possible model is to assume that \(\vect{b}\) should have constant entries, so the matrix \(A\) is a single columns of all ones. In this case, the normal equations provide \(\hat{\vect{x}}\) as the vector with a single entry that is the average of the entries of the vector \(\vect{b}\text{.}\) If we write the entries of \(\vect{b}\) as \(b_i\) and write the average as \(\overline{b}\text{,}\) then the length of the residual vector, squared, is

When we use a more sophisticated model, the residual vector will be shorter, and if we denote the entries of \(\hat{\vect{b}}\) by \(\hat{b}_i\text{,}\) we have

The notation for these expressions is read as a “sum of squares”, either “total” or “residual”. Then the coefficient of determination is defined as

When we actually have a solution to the original system, then the residual vector for the model is the zero vector, so \(\mathrm{SS}_\mathrm{res}=0\) and \(R^2=1\text{.}\) When the model does not do better than the simplest possible model, then \(\mathrm{SS}_\mathrm{res} = \mathrm{SS}_\mathrm{tot}\) and \(R^2=0\text{.}\) So the hope is that a model will fit the data with an \(R^2\) as close to \(1\) as possible. What is considered “close” will vary from situation to situation. What is considered close for economics data may not be considered close for data from a physics experiment.

Subsection 4.1.3 Computing Least-Squares Solutions

In Exercise Checkpoint 4.1.5 we suggest formulating general expressions for least-squares estimates of three parameters of a model. This begins to get a bit absurd when you invert a matrix larger than size \(2\text{.}\) Instead it makes sense to formulate from the data the matrix \(\adjoint{A}A\) and the vector \(\adjoint{A}\vect{b}\text{,}\) and then proceed to a numerical solution. As \(\adjoint{A}A\) is a Hermitian positive definite matrix (Theorem 1.8.2, it can be decomposed into a Cholesky factorization twice as fast as we can decompose an arbitrary matrix into an LU factorization [provisional cross-reference: cost to form Cholesky]. The Cholesky factorization allows a round of forward-solving followed by a round of back-solving to find a solution to the normal equations.

But it gets better. Consider solving the system \(A\vect{x}=\vect{b}\text{.}\) Suppose \(A\) has rank \(n\) and we have a thin QR decomposition of \(A\text{,}\) \(A=QR\) where \(Q\) has orthogonal columns and \(R\) is an upper triangular square matrix with positive diagonal entries. Notice in particular that \(R\) is invertible. Then

Since \(R\) is invertible, so is \(\adjoint{R}\text{,}\) and we can multiply both sides of the above equation by \(\inverse{\left(\adjoint{R}\right)}\) to obtain the system

This new system has an upper triangular coefficient matrix, so one round of back-solving will provide a solution for \(\hat{\vect{x}}\text{.}\) And the product \(\adjoint{Q}\vect{b}\) will be reasonably well-behaved due to the orthogonal columns of \(Q\text{.}\) Notice that in the case of a square matrix \(A\) (of full rank), the matrix \(Q\) will be invertible, and the same system that in general leads to the least-squares solution will also provide the exact solution in this special case.

This is an excellent example of how knowledge of a specialized situation, and a matrix decomposition, can be parlayed into better accuracty and greater speed. But it is not free, there is still significant effort required to compute the QR decomposition in the first place.