38.2: Finding the best solution in an overdetermined system

%matplotlib inline
import matplotlib.pylab as plt
import numpy as np
import sympy as sym
import time
sym.init_printing(use_unicode=True)

Let \(Ax=y\) be a system of \(m\) linear equations in \(n\) variables. A least squares solution of \(Ax=y\) is an solution \(\hat{x}\) in \(R^n\) such that:

\[ \min_{\hat{x}}\|y - A\hat{x}\| \nonumber \]

Note we substitute \(y\) for our typical variable \(b\) here because we will use \(b\) later to represent the intercept to a line and we want to try and avoid confusion in notation. It also consistent with the picture above.

In other words, \(\hat{x}\) is a value of \(x\) for which \(Ax\) is as close as possible to \(y\). From previous lectures, we know this to be true if the vector \(y-A\hat{x}\) is orthogonal (perpendicular) to the column space of \(A\).

We also know that the dot product is zero if two vectors are orthogonal. So we have \(a \cdot (Ax - y) = 0 \), for all vectors \(a\) in the column space of \(A\).

The columns of \(A\) span the column space of \(A\). Denote the columns of \(A\) as \(A = [a_1, \cdots, a_n]\). Then we have \(a_1^\top (Ax - y) = 0\), it is the same as taking the transpose of \(A\) and doing a matrix multiply \(A^\top (Ax - y) = 0\).

\[a_1^\top (Ax - y) = 0, \\ a_2^\top(Ax-y)=0\\\vdots \\a_n^\top(Ax-y)=0 \nonumber \]

That is:

\[A^\top Ax = A^\top y \nonumber \]

The above equation is called the least squares solution to the original equation \(Ax=y\). The matrix \(A^\top A\) is symmetric and invertable. Then solving for \(\hat{x}\) can be calculated as follows:

\[x = (A^\top A)^{-1}A^\top y \nonumber \]

The matrix \((A^\top A)^{-1}A^\top\) is also called the left inverse.

Example

A researcher has conducted experiments of a particular Hormone dosage in a population of rats. The table shows the number of fatalities at each dosage level tested. Determine the least squares line and use it to predict the number of rat fatalities at hormone dosage of 22.

Hormone level	20	25	30	35	40	45	50
Fatalities	101	115	92	64	60	50	49

H = [20,25,30,35,40,45,50]
f = [101,115, 92,64,60,50,49]

plt.scatter(H,f)
plt.xlabel('Hormone Level')
plt.ylabel('Fatalities')
f = np.matrix(f).T

We want to determine a line that is expressed by the following equation

\[f = aH + b, \nonumber \]

to approximate the connection between Hormone dosage (\(H\)) and Fatalities \(f\). That is, we want to find \(a\) (slope) and \(b\) (y-intercept) for this line. First we define the variable \(
x = \left[
\begin{matrix}
a \\
b
\end{matrix}
\right]
\) as the column vector that needs to be solved.

Do This

Rewrite the system of equations to the form \(Ax=y\) by defining your numpy matrices A and y using the data from above:

#put your code here

Question

Calculate the square matrix \(C = A^\top A\) and the modified right hand side vector as \(A^\top y\) (Call it Aty):

#put your code here

Question

Find the least squares solution by solving \(Cx=A^\top y\) for \(x\).

# Put your code here

Question

Given the solution above, define the two scalars slope a and y-intercept b.

#put your code here

The following code will Plot the original data and the line estimated by the coefficients found in the above quation.

H = [20,25,30,35,40,45,50]
f = [101,115, 92,64,60,50,49]
plt.scatter(H,f)

H2 = np.linspace(np.min(H), np.max(H))

f2 = a * H2 + b

plt.plot(H2, f2)

Question

Repeat the above analysis but now with a eight-order polynomial.

Question

Play with the interactive function below by adjusting the degree of the least-squares fit approximation. Then extend the x_min and x_max parameters. Do you think that an eight-order polynomial is a good model for this dataset? Why or why not?

from ipywidgets import interact, fixed
import ipywidgets as widgets

@interact(x=fixed(H), y=fixed(f), degree=widgets.IntSlider(min=1, max=8, step=1, value=8), x_min=widgets.IntSlider(min=min(H)-10, max=min(H), step=1, value=min(H)), x_max=widgets.IntSlider(min=max(H), max=max(H)+10, step=1, value=max(H)))
def graphPolyN(x, y, x_min, x_max, degree):
    p = np.polyfit(x, y, degree)
    f = np.poly1d(p)
    
    x_pred = np.linspace(x_min, x_max, 1000)
    y_pred = f(x_pred)
    
    plt.scatter(x, y, color="red")
    plt.plot(x_pred, y_pred)

Question

Check the rank of \(C=A^\top A\) for the previous case. What do you get? Why?