Lab 9 - Gradient Descent¶
This lab focuses on gradient descent and related concepts. You'll need the following libraries
## Standard Libraries
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import sklearn
import math
Part 1 - Squared Error Cost¶
For many structured machine learning models, the goal is to learn the values of weight parameters in the model's structure that will minimize a certain cost function. We'll begin by considering the squared error cost function: $$\text{Cost} = \frac{1}{n}(\mathbf{y} - \hat{\mathbf{y}})^T(\mathbf{y} - \hat{\mathbf{y}})$$
Here $\mathbf{y}$ is an $n \times 1$ vector of observed outcomes, and $\hat{\mathbf{y}}$ is an $n \times 1$ vector of predicted outcomes, which are found via a function involving $\mathbf{X}$, the $n \times p$ data matrix, and the model's weight parameters.
For linear regression, $\hat{\mathbf{y}} = \mathbf{X}\hat{\mathbf{w}}$, where $\hat{\mathbf{w}}$ is an $p \times 1$ vector of learned weights. Below we implement our own squared error cost function for linear regression:
## Squared Error Cost for Linear Regression
def cost_function(X, y, w):
n = len(y)
pred_y = np.dot(X,w)
err = y - pred_y
cost = (1./n)*np.dot(np.transpose(err),err)
return cost
The question below seeks to build intuition on how the learned weights of a model relate to the cost function for a simple simulated dataset created by the code below:
## Simulated Data for Question 1
np.random.seed(0)
x = np.linspace(-10, 10, 100)
y = 3 * x + np.random.normal(0, 3, 100)
Question #1:
- Part A: Consider the simulated data provided above. What are the dimensions of $\mathbf{X}$, $\mathbf{y}$, and $\hat{w}$ if a simple linear regression model with a single predictor and no intercept is used?
- Part B: Modify the provided cost function to account for $\mathbf{X}$ being a one-dimensional vector and $w_1$ being a scalar in the simulated dataset and corresponding model. Then create a sequence of 100 values for $\hat{w}_1$ ranging from 2 to 4, then display a line graph that shows the squared error cost associated with each value of $\hat{w}_1$.
- Part C: Suppose the current iteration of a gradient descent algorithm is using $\hat{w}_1 = 2$. Relying on your graph from Part B and without calculating anything, what is the sign of the gradient at this value of $\hat{w}_1$? How does this relate to the way that $\hat{w}_1$ should be adjusted to improve the cost?
- Part D: Now consider $\hat{w}_1 = 2.7$. Again relying on your graph from Part B, is the magnitude of the gradient smaller or larger than it was for $\hat{w}_1 = 2$ in Part C? Briefly explain how you arrived at your answer.
- Part E: Early in the semester we introduced the concept of irreducible error. What is the irreducible error in this application? How does it relate to the graph you created in Part B?
Part 2 - Gradient Descent¶
Gradient descent iteratively improves the parameters of a model (its weights) using the following update scheme:
$$ \hat{\mathbf{w}}^{(j)} = \hat{\mathbf{w}}^{(j-1)} - \alpha \cdot \frac{\partial \text{Cost}\!\left(\mathbf{w}^{(j-1)}\right)} {\partial \mathbf{w}} $$
For linear regression, the gradient (in matrix notation) is:
$$ \nabla \text{Cost} = \frac{-2}{n}\mathbf{X}^T \left( \mathbf{y} - \mathbf{X}\hat{\mathbf{w}} \right) $$
Note that you should be comfortable deriving this gradient yourself, either using substitution and chain rule, or by expanding the cost function.
Provided below is a simplistic gradient descent function that performs a fixed number of iterations:
## Define gradient descent function for linear regression
def grad_descent(X, y, w, alpha, n_iter):
costs = np.zeros(n_iter)
n = len(y)
for i in range(n_iter):
pred_y = np.dot(X,w)
grad = (-2./n)*np.dot(np.transpose(X), (y - pred_y))
w = w - alpha*grad
costs[i] = cost_function(X, y, w)
return w, costs
We can try this function out on the Iowa City Home Sales data set, using all of the numeric features to predict sale.amount:
## Run gradient descent (for the original IC homes data)
from sklearn.preprocessing import StandardScaler
ic = pd.read_csv("https://remiller1450.github.io/data/IowaCityHomeSales.csv")
ic_y = ic['sale.amount']
ic_X = ic.select_dtypes("number").drop('sale.amount',axis=1)
ic_Xs = StandardScaler().fit_transform(ic_X)
gdres = grad_descent(X=ic_Xs,y=ic_y,w=np.zeros(12),alpha=0.1, n_iter=100)
## Plot costs by iteration
plt.plot(np.linspace(0, 100, 100),gdres[1])
plt.xlabel("Epoch")
plt.ylabel("Cost")
plt.show()
We can see that the cost is leveling out towards the end of the plot (100th iteration of gradient descent). We can check how the weight estimates at this iteration compare to those obtained via the sklearn implementation of linear regression:
from sklearn.linear_model import LinearRegression
lmf = LinearRegression(fit_intercept=False).fit(ic_Xs, ic_y)
pd.DataFrame({"Feature": ic_X.columns, "sklearn Est": lmf.coef_, "Gradient Descent Est": gdres[0]})
| Feature | sklearn Est | Gradient Descent Est | |
|---|---|---|---|
| 0 | built | -1571.295331 | -2083.118805 |
| 1 | bedrooms | 2584.321728 | 1076.373337 |
| 2 | area.base | -225.923612 | -235.383234 |
| 3 | area.add | -5778.162652 | -6152.012101 |
| 4 | area.bsmt | -608.856265 | 1687.308121 |
| 5 | area.garage1 | -3022.003497 | -3088.592269 |
| 6 | area.garage2 | 1522.143087 | 1724.586773 |
| 7 | area.living | 159.080164 | 5188.784893 |
| 8 | area.lot | 1813.221752 | 1997.703559 |
| 9 | lon | 181.642943 | 310.978372 |
| 10 | lat | -520.807500 | -466.934918 |
| 11 | assessed | 90251.732347 | 85369.249327 |
Some of our weight estimates are close, but we aren't quite there yet. We can try again using more iterations of gradient descent:
## Try 2000 iterations of gradient descent
gdres = grad_descent(X=ic_Xs,y=ic_y,w=np.zeros(12),alpha=0.1, n_iter=2000)
pd.DataFrame({"Feature": ic_X.columns, "sklearn Est": lmf.coef_, "Gradient Descent Est": gdres[0]})
| Feature | sklearn Est | Gradient Descent Est | |
|---|---|---|---|
| 0 | built | -1571.295331 | -1571.295331 |
| 1 | bedrooms | 2584.321728 | 2584.321728 |
| 2 | area.base | -225.923612 | -225.923612 |
| 3 | area.add | -5778.162652 | -5778.162652 |
| 4 | area.bsmt | -608.856265 | -608.856265 |
| 5 | area.garage1 | -3022.003497 | -3022.003497 |
| 6 | area.garage2 | 1522.143087 | 1522.143087 |
| 7 | area.living | 159.080164 | 159.080164 |
| 8 | area.lot | 1813.221752 | 1813.221752 |
| 9 | lon | 181.642943 | 181.642943 |
| 10 | lat | -520.807500 | -520.807500 |
| 11 | assessed | 90251.732347 | 90251.732347 |
Here we obtained the least squares solution exactly. But you might be wondering if we could have simply used a larger learning rate instead of going up to 2,000 iterations. Below is the cost curve when the learning rate is set to $\alpha = 0.3$ and 100 iterations of gradient descent are performed:
## What happens with a larger learning rate?
gdres = grad_descent(X=ic_Xs,y=ic_y,w=np.zeros(12),alpha=0.3, n_iter=100)
plt.plot(np.linspace(0, 100, 100),gdres[1])
plt.xlabel("Epoch")
plt.ylabel("Cost")
plt.show()
With $\alpha = 0.3$ we are overshooting the minimum of the cost function by such a large degree that each new set of weights actually increases the cost, a problem that compounds as we perform more and more iterations.
Question #2: Consider a simple model involving a single predictor and single weight: $$\mathbf{y} = exp(w_1 \cdot \mathbf{x}) + \mathbf{\epsilon}$$
- Part A: Write out the squared-error cost function and differentiate with respect to $w_1$ using chain rule.
- Part B: Modify the functions provided earlier in this lab to perform gradient descent using the gradient you calculated in Part A.
- Part C: Apply your gradient descent function to the simulated data provided below, performing enough iterations to reach a stable estimate of $w_1$. Show your results using a line graph displaying the cost by epoch.
- Part D: Compare your final estimate of $w_1$ with the true value used to generate the data. Did gradient descent seem to effectively estimate this parameter?
## Simulated Data for Question 2
np.random.seed(0)
x = np.linspace(-1.5, 1.5, 100)
y = np.exp(-2 * x) + np.random.normal(0, 1, 100)
Part 3 - Mini-Batch Gradient Descent¶
Most modern machine learning models are trained using mini-batch gradient descent, where the training set is broken into batches and a gradient descent update is performed separately for each batch.
The code below modifies the gradient descent function provided earlier in the lab to perform mini-batch gradient descent:
def batch_grad_descent(X, y, w, alpha, n_iter, batch_size):
costs = np.zeros(n_iter)
n = len(y)
for epoch in range(n_iter):
# Shuffle observations
idx = np.random.permutation(n)
X_shuff = X[idx, :]
y_shuff = y.iloc[idx].values
# Process all the mini-batches
for start in range(0, n, batch_size):
end = min(start + batch_size, n)
batch_X = X_shuff[start:end, :]
batch_y = y_shuff[start:end]
pred_y = np.dot(batch_X, w)
grad = (-2.0 / len(batch_y)) * np.dot(batch_X.T,(batch_y - pred_y))
w = w - alpha * grad
# Cost on full dataset after the epoch
costs[epoch] = cost_function(X, y, w)
return w, costs
We can try out this approach on the Iowa City home sales example:
bgdres = batch_grad_descent(X=ic_Xs,y=ic_y,w=np.zeros(12),batch_size = 5, alpha=0.001, n_iter=500)
plt.plot(np.linspace(0, 500, 500),bgdres[1])
plt.xlabel("Epoch")
plt.ylabel("Cost")
plt.show()
Notice we get slightly noisy convergence to a similar minimum value of the cost function. We can further investigate by comparing the final weights (after 500 epochs) to the weights estimated by sklearn:
## Compare weights
pd.DataFrame({"Feature": ic_X.columns, "sklearn Est": lmf.coef_, "Batch Gradient Descent Est": bgdres[0]})
Question #3:
- Part A: In the Iowa City home sales example shown above the final weight estimates after 500 training epochs are somewhat different from the weights estimated by
sklearn. If we increased the number of training epochs, would we expect our weight estimates to exactly match those found bysklearn? Briefly explain. - Part B: Modify the functions you created in Question 2 to perform mini-batch gradient descent. Test your function using batch sizes of 4, 20, and 50, displaying the learning curve for each.
- Part C: Using the learning curves from Part B, which batch size produced the noisiest optimization path? Which appeared to converge most smoothly? Explain why these differences occur and describe one advantage and one disadvantage of using a smaller batch size. Hint: Because an initial value of 0 is so far from the optimal value of $w_1$, the learning curves in Part B tend to be dominated by the large costs at the start of the gradient descent. Consider changing the initial value to $-1.7$ to better see the differences across batch sizes.