6.5: Orthogonal Least Squares (2024)

Activity 6.5.2.

Suppose we have three data points \((1,1)\text{,}\) \((2,1)\text{,}\) and \((3,3)\) and that we would like to find a line passing through them.

1. Plot these three points in Figure 6.5.2. Are you able to draw a line that passes through all three points?

   

• Let's write the conditions that would describe a line passing through the points. Remember that the equation of a line can be written as \(b + mx=y\) where \(m\) is the slope and \(b\) is the \(y\)-intercept. We will try to find \(b\) and \(m\) so that the three points lie on the line.

  The first data point \((1,1)\) gives an equation for \(b\) and \(m\text{.}\) In particular, we know that when \(x=1\text{,}\) then \(y=1\) so we have \(b + m(1) = 1\) or \(b + m = 1\text{.}\) Use the other two data points to create a linear system describing \(m\) and \(b\text{.}\)

• We have obtained a linear system having three equations, one from each data point, for the two unknowns \(b\) and \(m\text{.}\) Identify a matrix \(A\) and vector \(\mathbf b\) so that the system has the form \(A\mathbf x=\mathbf b\text{,}\) where \(\mathbf x=\twovec bm\text{.}\)

  Notice that the unknown vector \(\mathbf x=\twovec bm\) describes the line that we seek.

• Is there a solution to this linear system? How does this question relate to your attempt to draw a line through the three points above?
• Since this system is inconsistent, we know that \(\mathbf b\) is not in the column space \(\col(A)\text{.}\) Find an orthogonal basis for \(\col(A)\) and use it to find the orthogonal projection \(\widehat\mathbf b\) of \(\mathbf b\) onto \(\col(A)\text{.}\)
• Since \(\widehat\mathbf b\) is in \(\col(A)\text{,}\) the equation \(A\mathbf x = \widehat\mathbf b\) is consistent. Find its solution \(\mathbf x = \twovec{b}{m}\) and sketch the line \(y=b + mx\) in Figure 6.5.2. We say that this is the line of best fit.

This activity illustrates the idea behind a technique known as orthogonal least squares, which we have been working toward throughout this chapter. If the data points are denoted as \((x_i, y_i)\text{,}\) we construct the matrix \(A\) and vector \(\mathbf b\) as

\begin{equation*} A = \begin{bmatrix} 1 & x_1 \\ 1 & x_2 \\ 1 & x_3 \\ \end{bmatrix},\hspace{24pt} \mathbf b = \threevec{y_1}{y_2}{y_3}\text{.} \end{equation*}

With the vector \(\mathbf x=\twovec bm\) representing the line \(b+mx = y\text{,}\) we see that the equation \(A\mathbf x=\mathbf b\) describes a line passing through all the data points. In our example, it is visually apparent that there is no such line, a fact confirmed by the inconsistency of the equation \(A\mathbf x=\mathbf b\text{.}\)

Remember that \(\bhat\text{,}\) the orthogonal projection of \(\mathbf b\) onto \(\col(A)\text{,}\) is the closest vector in \(\col(A)\) to \(\mathbf b\text{.}\) Therefore, when we solve the equation \(A\mathbf x=\bhat\text{,}\) we are finding the vector \(\mathbf x\) so that \(A\mathbf x = \threevec{b+mx_1}{b+mx_2}{b+mx_3}\) is as close to \(\mathbf b=\threevec{y_1}{y_2}{y_3}\) as possible. Let's think about what this means within the context of this problem.

The difference \(\mathbf b-A\mathbf x = \threevec{y_1-(b+mx_1)}{y_2-(b+mx_2)}{y_3-(b+mx_3)}\) so that the square of the distance between \(A\mathbf x\) and \(\mathbf b\) is

\begin{align*} \len{\mathbf b - A\mathbf x}^2 & =\\ & \left(y_1-(b+mx_1)\right)^2 + \left(y_2-(b+mx_2)\right)^2 + \left(y_3-(b+mx_3)\right)^2\text{.} \end{align*}

Our approach finds the values for \(b\) and \(m\) that make this sum of squares as small as possible, which explains why we call this a least squares problem.

Drawing the line defined by the vector \(\mathbf x=\twovec bm\text{,}\) the quantity \(y_i - (b + mx_i)\) reflects the vertical distance between the line and the data point \((x_i, y_i)\text{,}\) as shown in Figure 6.5.5. Seen in this way, the square of the distance \(\len{\mathbf b-A\mathbf x}^2\) is a measure of how much the line defined by the vector \(\mathbf x\) misses the data points. The solution to the least squares problem is the line that misses the data points by the smallest amount possible.

Activity 6.5.3.

The rate at which a cricket chirps is related to the outdoor temperature, as reflected in some experimental data that we'll study in this activity. The chirp rate \(C\) is expressed in chirps per second while the temperature \(T\) is in degrees Fahrenheit. Evaluate the following cell to load in the data:

Evaluating this cell also provides:

• the vectors chirps and temps formed from the columns of the dataset.
• the command onesvec(n), which creates an \(n\)-dimensional vector whose entries are all one.
• Remember that you can form a matrix whose columns are the vectors v1 and v2 with matrix([v1, v2]).T.

We would like to represent this relationship by a linear function

\begin{equation*} \beta_0 + \beta_1 C = T\text{.} \end{equation*}

1. Use the first data point \((C_1,T_1)=(20.0,88.6)\) to write an equation involving \(\beta_0\) and \(\beta_1\text{.}\)
2. Suppose that we represent the unknowns using a vector \(\mathbf x = \twovec{\beta_0}{\beta_1}\text{.}\) Use the 15 data points to create the matrix \(A\) and vector \(\mathbf b\) so that the linear system \(A\mathbf x= \mathbf b\) describes the unknown vector \(\mathbf x\text{.}\)
3. Write the normal equations \(A^TA\xhat = A^T\mathbf b\text{;}\) that is, find the matrix \(A^TA\) and the vector \(A^T\mathbf b\text{.}\)
4. Solve the normal equations to find \(\xhat\text{,}\) the least squares approximate solution to the equation \(A\mathbf x=\mathbf b\text{.}\) Call your solution xhat since x has another meaning in Sage.

   What are the values of \(\beta_0\) and \(\beta_1\) that you found?

5. If the chirp rate is 22 chirps per second, what is your prediction for the temperature?

   You can plot the data and your line, assuming you called the solution xhat, using the cell below.

This example demonstrates an approach called linear regression, in which a collection of data is modeled using a linear function found by solving a least squares problem. Once we have the linear function that best fits the data, we can make predictions about situations that we haven't encountered in the data.

If we're going to use our function to make predictions, it's natural to ask how much confidence we have in these predictions. This is a statistical question that leads to a rich and well-developed theory, which we won't explore in much detail here. However, there is one simple measure of how well our linear function fits the data that is known as the coefficient of determination and denoted by \(R^2\text{.}\)

We have seen that the square of the distance \(\len{\mathbf b-A\mathbf x}^2\) measures the amount by which the line fails to pass through the data points. When the line is close to the data points, we expect this number to be small. However, the size of this measure depends on the scale of the data. For instance, the two lines shown in Figure 6.5.8 seem to fit the data equally well, but \(|\mathbf b-A\xhat|^2\) is 100 times larger on the right.

The coefficient of determination \(R^2\) is defined by normalizing \(|\mathbf b-A\xhat|^2\) so that it is independent of the scale. Recall that we described how to demean a vector in Section 6.1: given a vector \(\mathbf v\text{,}\) we obtain \(\widetilde{\mathbf v}\) by subtracting the average of the components from each component.

A more complete explanation of this definition relies on the concept of variance, which we explore in Exercise 6.5.6.11 and the next chapter. For the time being, it's enough to know that \(0\leq R^2 \leq 1\) and that the closer \(R^2\) is to 1, the better the line fits the data. In our original example, illustrated in Figure 6.5.8, we find that \(R^2 = 0.75\text{,}\) and in our study of cricket chirp rates, we have \(R^2=0.69\text{.}\) However, assessing the confidence we have in predictions made by solving a least squares problem can require considerable thought, and it would be naive to rely only on the value of \(R^2\text{.}\)

Activity 6.5.4.

1. Suppose we are interested in finding the least squares approximate solution to the equation \(A\mathbf x = \mathbf b\) and that we have the \(QR\) factorization \(A=QR\text{.}\) Explain why the least squares approximation solution is given by solving

   \begin{align*} A\xhat & = QQ^T\mathbf b \\\\ QR\xhat & = QQ^T\mathbf b \\ \end{align*}

2. Multiply both sides of the second expression by \(Q^T\) and explain why

   \begin{equation*} R\xhat = Q^T\mathbf b. \end{equation*}

   Since \(R\) is upper triangular, this is a relatively simple equation to solve using back substitution, as we saw in Section 5.1. We will therefore write the least squares approximate solution as

   \begin{equation*} \xhat = R^{-1}Q^T\mathbf b, \end{equation*}

   and put this to use in the following context.

3. Brozak’s formula, which is used to calculate a person's body fat index \(BFI\text{,}\) is

   \begin{equation*} BFI = 100 \left(\frac{4.57}{\rho} - 4.142\right) \end{equation*}

   where \(\rho\) denotes a person's body density in grams per cubic centimeter. Obtaining an accurate measure of \(\rho\) is difficult, however, because it requires submerging the person in water and measuring the volume of water displaced. Instead, we will gather several other body measurements, which are more easily obtained, and use it to predict \(BFI\text{.}\)

   For instance, suppose we take 10 patients and measure their weight \(w\) in pounds, height \(h\) in inches, abdomen \(a\) in centimeters, wrist circumference \(r\) in centimeters, neck circumference \(n\) in centimeters, and \(BFI\text{.}\) Evaluating the following cell loads and displays the data.

   In addition, that cell provides:
   1. vectors weight, height, abdomen, wrist, neck, and BFI formed from the columns of the dataset.
   2. the command onesvec(n), which returns an \(n\)-dimensional vector whose entries are all one.
   3. the command QR(A) that returns the \(QR\) factorization of \(A\) as Q, R = QR(A).
   4. the command demean(v), which returns the demeaned vector \(\widetilde{\mathbf v}\text{.}\)

   We would like to find the linear function

   \begin{equation*} \beta_0 + \beta_1w + \beta_2h + \beta_3a + \beta_4r + \beta_5n = BFI \end{equation*}

   that best fits the data.

   Use the first data point to write an equation for the parameters \(\beta_0,\beta_1,\ldots,\beta_5\text{.}\)

4. Describe the linear system \(A\mathbf x = \mathbf b\) for these parameters. More specifically, describe how the matrix \(A\) and the vector \(\mathbf b\) are formed.
5. Construct the matrix \(A\) and find its \(QR\) factorization in the cell below.
6. Find the least squares approximate solution \(\xhat\) by solving the equation \(R\xhat = Q^T\mathbf b\text{.}\) You may want to use N(xhat) to display a decimal approximation of the vector. What are the parameters \(\beta_0,\beta_1,\ldots,\beta_5\) that best fit the data?
7. Find the coefficient of determination \(R^2\) for your parameters. What does this imply about the quality of the fit?
8. Suppose a person's measurements are: weight 190, height 70, abdomen 90, wrist 18, and neck 35. Estimate this person's \(BFI\text{.}\)

To summarize, we have seen that

Activity 6.5.5.

1. Suppose that we have a small dataset containing the points \((0,2)\text{,}\) \((1,1)\text{,}\) \((2,3)\text{,}\) and \((3,3)\text{,}\) such as appear when the following cell is evaluated. In addition to loading and plotting the data, evaluating that cell provides the following commands:
   • Q, R = QR(A) returns the \(QR\) factorization of \(A\text{.}\)
   • demean(v) returns the demeaned vector \(\widetilde{\mathbf v}\text{.}\)

   Let's fit a quadratic function of the form

   \begin{equation*} \beta_0 + \beta_1 x + \beta_2 x^2 = y \end{equation*}

   to this dataset.

   Write four equations, one for each data point, that describe the coefficients \(\beta_0\text{,}\) \(\beta_1\text{,}\) and \(\beta_2\text{.}\)

2. Express these four equations as a linear system \(A\mathbf x = \mathbf b\) where \(\mathbf x = \threevec{\beta_0}{\beta_1}{\beta_2}\text{.}\)

   Find the \(QR\) factorization of \(A\) and use it to find the least squares approximate solution \(\xhat\text{.}\)

3. Use the parameters \(\beta_0\text{,}\) \(\beta_1\text{,}\) and \(\beta_2\) that you found to write the quadratic function that fits the data. You can plot this function, along with the data, by entering your function in the appropriate place below.
4. What is your predicted \(y\) value when \(x=1.5\text{.}\)
5. Find the coefficient of determination \(R^2\) for the quadratic function? What does this say about the quality of the fit?
6. Now fit a cubic polynomial of the form

   \begin{equation*} \beta_0 + \beta_1x + \beta_2 x^2 + \beta_3x^3 = y \end{equation*}

   to this dataset.

7. Find the coefficient of determination \(R^2\) for the cubic function. What does this say about the quality of the fit?
8. What do you notice when you plot the cubic function along with the data? How does this reflect the value of \(R^2\) that you found?

The matrices \(A\) that you created in the last activity when fitting a quadratic and cubic function to a dataset have a special form. In particular, if the data points are labeled \((x_i, y_i)\) and we seek a degree \(k\) polynomial, then

\begin{equation*} A = \begin{bmatrix} 1 & x_1 & x_1^2 & \ldots & x_1^k \\ 1 & x_2 & x_2^2 & \ldots & x_2^k \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_m & x_m^2 & \ldots & x_m^k \\ \end{bmatrix}. \end{equation*}

This is called a Vandermonde matrix of degree \(k\text{.}\)

Activity 6.5.6.

This activity explores a dataset describing Arctic sea ice and that comes from Sustainability Math.

Evaluating the cell below will plot the extent of Arctic sea ice, in millions of square kilometers, during the twelve months of 2012.

In addition, you have access to a few special variables and commands:

• month is the vector of month values and ice is the vector of sea ice values from the table above.
• vandermonde(x, k) constructs the Vandermonde matrix of degree \(k\) using the points in the vector x.
• Q, R = QR(A) provides the \(QR\) factorization of \(A\text{.}\)
• demean(v) returns the demeaned vector \(\widetilde{\mathbf v}\text{.}\)

1. Find the vector \(\xhat\text{,}\) the least squares approximate solution to the linear system that results from fitting a degree 5 polynomial to the data.
2. If your result is stored in the variable xhat, you may plot the polynomial and the data together using the following cell.
3. Find the coefficient of determination \(R^2\) for this polynomial fit.
4. Repeat these steps to fit a degree 8 polynomial to the data, plot the polynomial with the data, and find \(R^2\text{.}\)
5. Repeat one more time by fitting a degree 11 polynomial to the data, plotting it, and finding \(R^2\text{.}\)

   It's certainly true that higher degree polynomials fit the data better, as seen by the increasing values of \(R^2\text{,}\) but that's not always a good thing. For instance, when \(k=11\text{,}\) you may notice that the graph of the polynomial wiggles a little more than we would expect. In this case, the polynomial is trying too hard to fit the data, which usually contains some uncertainty, especially if it's obtained from measurements. The error built in to the data is called noise, and its presence means that we shouldn't expect our polynomial to fit the data perfectly. When we choose a polynomial whose degree is too high, we give the noise too much weight in the model, which leads to some undesirable behavior, like the wiggles in the graph.

   Fitting the data with a polynomial whose degree is too high is called overfitting, a phenomenon that can appear in many machine learning applications. Generally speaking, we would like to choose \(k\) large enough to capture the essential features of the data but not so large that we overfit and build the noise into the model. There are ways to determine the optimal value of \(k\text{,}\) but we won't pursue that here.

6. Choosing a reasonable value of \(k\text{,}\) estimate the extent of Arctic sea ice at month 6.5, roughly at the Summer Solstice.