<jupyter_start><jupyter_text>CS4618: Artificial Intelligence IThe Normal Equation for OLS Regression    Derek Bridge    School of Computer Science and Information Technology    University College Cork Initialization$\newcommand{\Set}[1]{\{1\}}$ $\newcommand{\Tuple}[1]{\langle1\rangle}$ $\newcommand{\v}[1]{\pmb{1}}$ $\newcommand{\cv}[1]{\begin{bmatrix}1\end{bmatrix}}$ $\newcommand{\rv}[1]{[1]}$ $\DeclareMathOperator{\argmax}{arg\,max}$ $\DeclareMathOperator{\argmin}{arg\,min}$ $\DeclareMathOperator{\dist}{dist}$$\DeclareMathOperator{\abs}{abs}$<jupyter_code>%load_ext autoreload
%autoreload 2
%matplotlib inline
import pandas as pd
import numpy as np

from ipywidgets import interactive

import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec

from sklearn.preprocessing import add_dummy_feature

import numpy.linalg as npla<jupyter_output><empty_output><jupyter_text>Linear Regression    How does linear regression find the best model?    Recap:                    The learner is given a labeled training set ($m \times n$ matrix $\v{X}$ and                 $m$-dimensional vector $\v{y}$) and it inserts an extra element $\v{x}_0^{(i)}$ in each                row $i$, all of which will be 1                $$\v{X} = \begin{bmatrix}                    \v{x}_0^{(1)} = 1 & \v{x}_1^{(1)} & \v{x}_2^{(1)} & \ldots & \v{x}_n^{(1)} \\                    \v{x}_0^{(2)} = 1 & \v{x}_1^{(2)} & \v{x}_2^{(2)} & \ldots & \v{x}_n^{(2)} \\                    \vdots            & \vdots        & \vdots        & \vdots & \vdots \\                    \v{x}_0^{(m)} = 1 & \v{x}_1^{(m)} & \v{x}_2^{(m)} & \ldots & \v{x}_n^{(m)} \\                          \end{bmatrix}\,\,\,\,                   \v{y} = \cv{y^{(1)}\\ y^{(2)}\\ \vdots\\ y^{(m)}}                 $$                          An $n$-dimensional column vector                  $\v{\beta} = \cv{\v{\beta}_0\\ \v{\beta}_1\\ \vdots\\ \v{\beta}_n}$                 gives us a hypothesis $h_{\v{\beta}}$ where                 $$\hat{y} = h_{\v{\beta}}(\v{x}) = \v{x}\v{\beta}$$             In the case of Ordinary Least Squares (OLS) regression, we must find the vector $\v{\beta}$                 that minimizes the following loss function                 $$J(\v{X}, \v{y}, \v{\beta}) = \frac{1}{2m}\sum_{i=1}^m(h_{\v{\beta}}(\v{x}^{i)}) - \v{y}^{(i)})^2$$                          In Machine Learning, the variables whose values we are trying to find (in this case, $\v{\beta}$)                are called the parameters.<jupyter_code># Loss function for OLS regression (assumes X contains all 1s in its first column)
def J(X, y, beta):
    return np.mean((X.dot(beta) - y) ** 2) / 2.0<jupyter_output><empty_output><jupyter_text>2D Visualization of $J$    Let's visualize $J$ using the Cork Property Prices Dataset as the training set.    For a 2D visualization, we'll assume that $\v{\beta}_j = 0$ for all $j$ except $j = 1$.                    In other words, we are pretending that floor area is the only relevant feature.                            Then we can plot $J$ on the vertical axis against different values of        $\v{\beta}_{1}$ on the horizontal axis.<jupyter_code># Use pandas to read the CSV file
df = pd.read_csv("../datasets/dataset_corkA.csv")

# Get the feature-values (just flarea) and the target values 
X = df[["flarea"]].values
y = df["price"].values
def show_linear_model_and_loss(beta1):
    fig = plt.figure(figsize=(8, 4)) 
    gs = gridspec.GridSpec(1, 2, width_ratios=[2, 1]) 
    
    ax0 = plt.subplot(gs[0])
    plt.title("Training set")
    plt.xlabel("Floor area (sq metres)")
    plt.xlim(-500, 500)
    plt.ylabel("Price (000 euros)")
    plt.ylim(0, 1000)

    ax1 = plt.subplot(gs[1])
    plt.title("Loss function")
    plt.xlabel("$beta1$")
    plt.xlim(-10, 10)
    plt.ylabel("$J$")
    plt.ylim(0, 1500000)
    
    fig.tight_layout()
    
    ax0.scatter(df["flarea"], df["price"], color = 'green')
    xvals = np.linspace(-500, 500, 3)
    ax0.plot(xvals, xvals * beta1, color = "blue")
    ax1.scatter(beta1, J(X, y, beta1), color = "red")

    plt.show()
    
interactive_plot = interactive(show_linear_model_and_loss, beta1=(-10,10,.1))
interactive_plot<jupyter_output><empty_output><jupyter_text>Another 2D Visualization of $J$    Instead of making manual adjustments, let's use a loop to try several values for $\v{\beta}_1$<jupyter_code>fig = plt.figure(figsize=(8, 4)) 
gs = gridspec.GridSpec(1, 2, width_ratios=[2, 1]) 
ax0 = plt.subplot(gs[0])
plt.title("Training set")
plt.xlabel("Floor area (sq metres)")
plt.xlim(-500, 500)
plt.ylabel("Price (000 euros)")
plt.ylim(0, 1000)
ax1 = plt.subplot(gs[1])
plt.title("Loss function")
plt.xlabel("$beta1$")
plt.xlim(-10, 10)
plt.ylabel("$J$")
plt.ylim(0, 1500000)
fig.tight_layout()
ax0.scatter(df["flarea"], df["price"], color = 'green')

xvals = np.linspace(-500, 500, 3)
beta1s = np.linspace(-10, 10, 21)

for beta1 in beta1s:
    ax0.plot(xvals, xvals * beta1, color = "blue")

ax1.scatter(beta1s, [J(X, y, beta1) for beta1 in beta1s], color = "red")

plt.show()<jupyter_output><empty_output><jupyter_text>The loss function is convex.    Informally, this means:                                    in 2D, it is u-shaped;                                        it has a unique minimum.                            This is no accident: it follows from the way the loss function has been defined. 3D Visualization of $J$    Let's visualize $J$ again using the Cork Property Prices Dataset as the training set.    This time, we'll assume that $\v{\beta}_j = 0$ for all $j$ except $j = 2$ and $j = 3$.                    In other words, we are pretending that the number of bedrooms and bathrooms are the only relevant                features.                            This will be a 3D plot with $J$ on the vertical axis against different values of $\v{\beta}_2$         and $\v{\beta}_3$ on the horizontal axes.<jupyter_code>X = df[["bdrms", "bthrms"]].values
y = df["price"].values 
fig = plt.figure() 
ax = fig.add_subplot(111, projection='3d')
ax.set_title("Loss function")
ax.set_xlabel("$beta_2$")
ax.set_ylabel("$beta_3$")
ax.set_zlabel("$J$")
xvals = np.linspace(-100, 200, 301)
yvals = np.linspace(-100, 200, 301)
xxvals, yyvals = np.meshgrid(xvals, yvals)
zs = np.array([J(X, y, [beta2, beta3]) for beta2, beta3 in zip(xxvals.flatten(), yyvals.flatten())])
zvals = zs.reshape(xxvals.shape)
ax.plot_surface(xxvals, yyvals, zvals)
plt.show()<jupyter_output><empty_output><jupyter_text>Another Visualization of $J$    Here is the same data that we had in 3D but on a different kind of plot, called a contour plot.    In effect, it flattens the diagram above.    The lines connect points that have the same value for $J$.<jupyter_code>fig = plt.figure(figsize = (6, 6)) 
plt.title("Loss function")
plt.xlabel("$beta_2$")
plt.ylabel("$beta_3$")
xvals = np.linspace(-100, 200, 301)
yvals = np.linspace(-100, 200, 301)
xxvals, yyvals = np.meshgrid(xvals, yvals)
zs = np.array([J(X, y, [beta2, beta3]) for beta2, beta3 in zip(xxvals.flatten(), yyvals.flatten())])
zvals = zs.reshape(xxvals.shape)
C = plt.contour(xxvals, yyvals, zvals, 15, colors = "black")
plt.clabel(C, inline=1, fontsize=10)
plt.show()<jupyter_output><empty_output><jupyter_text>The 3D visualization and the contour plot show that this too is convex.    Informally, this means:                                    in 3D it is bowl-shaped;                                        again there is a unique minimum.                         Ordinary Least Squares Linear Regression    How do we find the best values for $\v{\beta}$?    There are two methods:                    The Normal Equation; and            Gradient Descent.             The Normal Equation    The normal equation solves for $\v{\beta}$:        $$\v{\beta} = (\v{X}^T\v{X})^{-1}\v{X}^T\v{y}$$                    In other words, the normal equation gives us the                parameters that minimize the loss function.                            Where does it come from?                    Take the gradient of the loss function: $\frac{1}{m}\v{X}^T(\v{X}\v{\beta} - \v{y})$ (see next slide).            Set it to zero: $\frac{1}{m}\v{X}^T(\v{X}\v{\beta} - \v{y}) = \v{0}$ (in fact, a                $(n+1)$-dimensional vector of zeros).                        Then do some algebraic manipulation to get $\v{\beta}$ on the left-hand side:                $\v{\beta} = (\v{X}^T\v{X})^{-1}\v{X}^T\v{y}$                         Partial Derivatives    We need the gradient of the loss function with regards to each $\v{\beta}_j$.                    In other words, how much the loss will change if we change $\v{\beta}_j$ a little.            With respect to a particular $\v{\beta}_j$, it is called the partial derivative.                Without doing the calculus, the partial derivatives of $J(\v{X}, \v{y}, \v{\beta})$ with respect to        $\v{\beta}_j$ are        $$\frac{\partial J(\v{X}, \v{y}, \v{\beta})}{\partial\v{\beta}_j} =           \frac{1}{m}\sum_{i=1}^m(\v{x}^{(i)}\v{\beta}_j - \v{y}^{(i)}) \times \v{x}_j^{(i)}$$        The gradient vector, $\nabla_{\v{\beta}}J(\v{X}, \v{y}, \v{\beta})$, is a vector        of each partial derivative:        $$\nabla_{\v{\beta}}J(\v{X}, \v{y}, \v{\beta}) =             \cv{\frac{\partial J(\v{X}, \v{y}, \v{\beta})}{\partial\v{\beta}_0}\\               \frac{\partial J(\v{X}, \v{y}, \v{\beta})}{\partial\v{\beta}_1}\\               \vdots\\               \frac{\partial J(\v{X}, \v{y}, \v{\beta})}{\partial\v{\beta}_n}} =              \cv{\frac{1}{m}\sum_{i=1}^m(\v{x}^{(i)}\v{\beta}_0 - \v{y}^{(i)}) \times \v{x}_0^{(i)}\\                \frac{1}{m}\sum_{i=1}^m(\v{x}^{(i)}\v{\beta}_1 - \v{y}^{(i)}) \times \v{x}_1^{(i)}\\                \vdots\\                \frac{1}{m}\sum_{i=1}^m(\v{x}^{(i)}\v{\beta}_n - \v{y}^{(i)}) \times \v{x}_n^{(i)}}        $$          And there is a vectorized way to compute it:         $$\nabla_{\v{\beta}}J(\v{X}, \v{y}, \v{\beta}) = \frac{1}{m}\v{X}^T(\v{X}\v{\beta} - \v{y})$$      The Normal Equation in scikit-learn    The fit method of scikit-learn's LinearRegression class does        what we have described:                    It inserts the extra column of 1s.            It calculates $\v{\beta}$ using the Normal Equation.             The Normal Equation in numpy    For the hell of it, let's see how to implement it ourselves.    Again for the purposes of this explanation, we will use the entire dataset as our training set.                    We will learn later that using all the data for training is usually not                the right thing to do.<jupyter_code># Get the feature-values and the target values 
X = df[["flarea", "bdrms", "bthrms"]].values
y = df["price"].values

# Add the extra column to X
X = add_dummy_feature(X)

# Solve the normal equation
beta = npla.inv(X.T.dot(X)).dot(X.T).dot(y)

# Display beta
beta<jupyter_output><empty_output><jupyter_text>Essentially, one line of code!    But there's a problem:                    The normal equation requires that $\v{X}^T\v{X}$ has an inverse.            But it might not (see earlier lecture about matrices).                When we discussed this before, we mentioned that, in some cases, we can use the         pseudo-inverse instead.                    This is one of those cases.                So the more robust way of writing this program is:<jupyter_code># Solve the normal equation - but using the pseudo-inverse
beta = npla.pinv(X.T.dot(X)).dot(X.T).dot(y)

# Display beta
beta<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1648103028193526, 'fraction_numerical': 0.023146536721197354, 'mean_word_length': 2.587078651685393, 'pattern_counts': {'":': 0, '<': 33, '<?xml version=': 0, '>': 33, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '229328', 'license': 'MIT', 'n_tokens_mistral': 3992, 'n_tokens_neox': 3732, 'n_words': 1117, 'path': 'ai1/lectures/AI1_06_normal_equation.ipynb', 'repo_name': 'derekgbridge/artificial_intelligence', 'size': '11492'}