Least-Squares Approximation

Often an exact solution to a problem in applied mathematics is difficult or impossible to obtain. However, it is usually just as useful to find an approximation to a solution. In particular, finding “linear approximations" is a powerful technique in applied mathematics. One basic case is the situation where a system of linear equations has no solution, and it is desirable to find a “best approximation" to a solution to the system.

We begin by defining the “best approximation” in a natural way, and showing that computing the best approximation reduces to solving a related system of linear equations called the normal equations. Next, we demonstrate a common application where a collection of data points is approximated by a line (or a curve). We conclude this section by showing that \(QR\)-factorization provides us with a more efficient way to solve the normal equations and compute the best approximation.

Best Approximate Solutions

Let
\[A=\begin{bmatrix}3 & 1\\1 & 2\\1 & 2\end{bmatrix}\quad \text {and}\quad \vec {b}=\begin{bmatrix}2\\1\\3\end{bmatrix}\]
Consider the matrix equation \(A\vec {x}=\vec {b}\). A quick examination of the last two rows should convince you that this equation has no solutions. In other words, \(\vec {b}\) is not in the span of the columns of \(A\).

If \(\vec {z}\) were an exact solution to \(A\vec {x}=\vec {b}\), then \(\vec {b}-A\vec {z}\) would be \(\vec {0}\). Since the equation does not have a solution, we will attempt to find the next best thing to a solution by finding \(\vec {z}\) such that \(\norm {\vec {b}-A\vec {z}}\) is as small as possible. The quantity \(\norm {\vec {b}-A\vec {z}}\) is called the error.

The following GeoGebra interactive will help you understand the geometry behind finding \(\vec {z}\).

RIGHT-CLICK and DRAG to rotate the image for a better view.

Record your best guess for \(\vec {z}\) – you will have a chance to check your answer in Example 2.

What did you discover about the geometry of minimizing \(\norm {\vec {b}-A\vec {z}}\)? Select all that apply.

\(\vec {z}\) is orthogonal to the plane spanned by the columns of \(A\). \(\norm {\vec {b}-A\vec {z}}\) is orthogonal to \(\text {col}(A)\). \(\vec {b}-A\vec {z}\) is orthogonal to \(\text {col}(A)\). \(A\vec {z}\) is orthogonal to \(\text {col}(A)\). \(A\vec {z}\) is an orthogonal projection of \(\vec {b}\) onto \(\text {col}(A)\).

Our geometric observations will help us develop a method for finding \(\vec {z}\) .

Suppose \(A\) is an \(m\times n\) matrix, and \(\vec {b}\) is a column vector in \(\RR ^m\). Consider the matrix equation \(A\vec {x}=\vec {b}\). If this equation does not have a solution, we can attempt to find a best approximation by finding \(\vec {z}\) which minimizes the error, \(\norm {\vec {b}-A\vec {z}}\). The expression \(\norm {\vec {b}-A\vec {z}}\) is also sometimes called the residual. The error (or the residual) is given in terms of a vector norm. Recall that our definition of the norm involves the sum of squares of the vector components. When we minimize the norm, we minimize the sum of squares. This is why the method we are describing is often referred to as least squares. We will explore this idea further later in this section.

In the case when \(\text {col}(A)\) is a subspace of \(\RR ^3\), we can see geometrically that \(\vec {z}\) is the best approximation if and only if \(A\vec {z}\) is an orthogonal projection of \(\vec {b}\) onto \(\text {col}(A)\), and the error is the magnitude of \(\vec {b}-A\vec {z}\), as shown below.

[Picture]

What we observed above, holds in general. We will use this fact to find \(\vec {z}\).

Every vector in \(\text {col}(A)\) can be written in the form \(A\vec {x}\) for some \(\vec {x}\) in \(\RR ^m\). Our goal is to find \(\vec {z}\) such that \(A\vec {z}\) is the orthogonal projection of \(\vec {b}\) onto \(\text {col}(A)\). By Corollary ??, every vector \(A\vec {x}\) in \(\text {col}(A)\) is orthogonal to \(\vec {b}-A\vec {z}\).

Therefore, we have

\begin{eqnarray} A^T(\vec {b}-A\vec {z})&=&\vec {0}\nonumber \\ A^T\vec {b}-A^TA\vec {z}&=&\vec {0}\nonumber \\ A^TA\vec {z}&=&A^T\vec {b}\label {eq:normalForZ} \end{eqnarray}

Since \(\vec {b}-A\vec {z}\) is normal to the subspace \(\text {col}(A)\), we call the system of linear equations in (??) the normal equations for \(\vec {z}\). If \(A^TA\) is invertible, then we can write

\begin{equation}\label {eq:leastSquaresZ} \vec {z}=(A^TA)^{-1}A^T\vec {b} \end{equation}

We will return to the question of invertibility of \(A^TA\) in Theorem 4. For now, let’s revisit the problem posed in Exploration .

We now come back to the question of when \(A^TA\) is invertible.

Let \(A\) be a matrix with linearly independent columns. We will show that \(\left (A^TA\right )\vec {x}=\vec {0}\) has only the trivial solution. For \(\vec {x}\), a solution of \(A^TA\vec {x}=\vec {0}\), we have

\begin{eqnarray*} \norm {A\vec {x}}^2&=&(A\vec {x})\dotp (A\vec {x})\\ &=&\left (A\vec {x}\right )^TA\vec {x}\\ &=&\vec {x}^TA^TA\vec {x}\\ &=&\vec {x}^T\dotp \vec {0}=0 \end{eqnarray*}

Therefore \(A\vec {x}=\vec {0}\). By linear independence of the columns of \(A\) we conclude that \(\vec {x}=\vec {0}\).

We summarize our findings in the following theorem.

Application of Least Squares to Curve Fitting

In Curve Fitting, we discussed how to fit a function to a set of data points, so that the graph of the function passes through each of the points. We also discussed why doing so is sometimes impossible (two points lie one above the other), and may not even be desirable (overfitting). In this section we will learn how to approximate a collection of data points with a line (or a curve) that fits the “trend" of the points. We will start with data that fit a linear pattern.

Consider the points \((1,1)\), \((2, 3)\) and \((4,4)\). These points do not lie on a straight line, but they have a general upward linear trend. (Typically there would be many more points to consider, but we will limit our exploration to what we can do by hand.) Our goal is to find a line that fits these points as closely as possible.

We are looking for a function \(f\) of the form \(f(x)=ax+b\) such that the following infeasible system is satisfied as closely as possible

\[\begin{matrix}a(1)&+&b&=&1\\a(2)&+&b&=&3\\a(4)&+&b&=&4\end{matrix}\]

From the first part of this section we know how to find a best approximation. By Theorem 6, we have

\begin{eqnarray*} \vec {z}=\begin{bmatrix}a\\b\end{bmatrix}&=&\left (A^TA\right )^{-1}A^T\vec {b}\\ &=&\left (\begin{bmatrix}1&2&4\\1&1&1\end{bmatrix}\begin{bmatrix} 1&1\\2&1\\4&1\end{bmatrix}\right )^{-1}\begin{bmatrix}1&2&4\\1&1&1\end{bmatrix}\begin{bmatrix} 1\\3\\4 \end{bmatrix}\\ &=&\begin{bmatrix}13/14\\1/2\end{bmatrix} \end{eqnarray*}

According to our computations, the line that best fits the data is given by

\[f(x)=\frac {13}{14}x+\frac {1}{2}\]
Let’s take a look.

[Picture]

We found this fit by minimizing \(\norm {\vec {b}-A\vec {z}}\). We will now investigate the meaning of this expression in relation to the line and the data points.

\begin{equation} \vec {b}-A\vec {z}=\begin{bmatrix}1\\3\\4\end{bmatrix}-\begin{bmatrix}(13/14)(1)+1/2\\(13/14)(2)+1/2\\{(13/14)(4)+1/2}\end{bmatrix}\approx \begin{bmatrix} -0.43\\0.64\\-0.21 \end{bmatrix} \end{equation}
Observe that each entry of \(\vec {b}-A\vec {z}\) is the signed vertical distance between a particular point and the line.

Instead of computing the error, \(\norm {\vec {b}-A\vec {z}}\), we will compute \(\norm {\vec {b}-A\vec {z}}^2\) to avoid the square root.

\begin{equation} \norm {\vec {b}-A\vec {z}}^2=(-0.43)^2+0.64^2+(-0.21)^2\approx 0.64 \end{equation}

Minimizing \(\norm {\vec {b}-A\vec {z}}\) also minimizes \(\norm {\vec {b}-A\vec {z}}^2\). Therefore, what we have minimized is the sum of squares of the vertical distances between the data points and the line. The following GeoGebra interactive will help you explore this idea.

In Exploration we discovered that \(\norm {\vec {b}-A\vec {z}}^2\) is the sum of squares of vertical distances between the given data points and the proposed line. By minimizing \(\norm {\vec {b}-A\vec {z}}\), we minimize the sum of squares of vertical distances. This observation holds in general. Given a collection of points \((x_1, y_1), (x_2, y_2),\dots ,(x_n, y_n)\), finding a linear function of the form \(f(x)=ax+b\) that best fits the points we would find a best solution to the system

\[\begin{bmatrix}x_1&1\\x_2&1\\\vdots& \vdots \\x_n&1\end{bmatrix}\begin{bmatrix}a\\b\end{bmatrix}=\begin{bmatrix}y_1\\y_2\\\vdots \\y_n\end{bmatrix}\]

by minimizing

\[\norm {\begin{bmatrix}y_1-(ax_1+b)\\y_2-(ax_2+b)\\\vdots \\y_n-(ax_n+b)\end{bmatrix}}^2=(y_1-(ax_1+b))^2+(y_2-(ax_2+b))^2+\dots +(y_n-(ax_n+b))^2\]

A geometric interpretation of \(y_i-(ax_i+b)\) is shown below.

[Picture]

The line we obtain in this fashion is called a line of best fit or a trendline, and the method we used is referred to as the method of least squares.

We can apply the method of least squares to find best fitting non-linear functions.

For more information about doing least squares with Octave, please see Octave for Chapter 9.

Practice Problems

Find the best approximation to a solution to the system of equations.
\[\begin{matrix}3x&+&y&+&z&=&6\\2x&+&3y&-&z&=&1\\2x&-&y&+&z&=&0\\3x&-&3y&+&3z&=&8\end{matrix}\]
Enter answers in fraction form below.
\[x=\answer {\frac {-20}{12}},\quad y=\answer {\frac {46}{12}},\quad z=\answer {\frac {95}{12}}\]
Find a linear function of best fit for the given set of data points. Examine how well your line fits the points by typing the equation of the line into the Desmos window.
\[(2,4), (4,3), (7,2), (8,1)\]
Enter your answers in fraction form.
\[f(x)=\answer {\frac {-6}{13}}x+\answer {\frac {64}{13}}\]

Find a linear function of best fit for the given set of data points. Examine how well your line fits the points by typing the equation of the line into the Desmos window.
\[(-2, 3), (-1,1), (0,0), (1, -2), (2, -4)\]
\[f(x)=\answer {-1.7}x+\answer {-0.4}\]

Use Octave to find the least squares approximating quadratic function for the following data points.
\[(-2,1),(0,0),(3,2),(4,3)\]
Round your answers to three decimal places.
\[f(x)=\answer {0.194}x^2+\answer {-0.024}x+\answer {0.127}\]
For more information about doing least squares with Octave, please see Octave for Chapter 9.

Text Source

A portion of this section has been adapted from W. Keith Nicholson, Linear Algebra with Applications, Lyryx 2021-A, Open Edition, pp. 308-319.