Gradient Descent

Gradient Descent (GD) is a generic optimization algorithm that can be applied to a wide range of problems. Note that GD is not a model: it is a training algorithm — a way to find a model's parameters. In this lecture, we use it to train linear regression, the same model we trained with the normal equation in the Linear Regression lecture.

The basic idea of GD is to iteratively tweak model parameters to minimize the cost function (error function or loss function) such as MSE we learned before.

Imagine walking down a mountain in thick fog. You cannot see the valley, but you can feel which way the ground slopes under your feet. If you keep stepping in the steepest downhill direction, you eventually reach the bottom. Gradient descent does the same on the cost function.

GD measures the local gradient of the cost function with regard to the parameter vector $\mathbf{w}$ and goes in the direction of descending gradient.

Once the gradient is zero, the minimum is reached.

The basic idea:

• random initialization: fill $\mathbf{w}$ and $b$ with random values
• take one small step (measured by the learning rate eta $\eta$) at a time trying to decrease the cost function (MSE)
• repeat till reaching the minimum (algorithm convergence) or approach arbitrarily close to the minimum

As shown in the figure below:

The code for the figures on this page is in the Appendix.

Learning Rate $\eta$ (eta)

If we set the learning rate too small, it may take a long time to converge:

If we set the learning rate too large, the algorithm may jump around and even diverge:

When the cost function is not convex (not like a bowl, with holes, ridges, plateaus, etc.), the algorithm may find only the local minimum or take a long time to cross the plateaus - we do a few random initialization hoping to land to the place that can lead us to the global minimum quickly.

NOTE: The cost function MSE is a convex function, which guarantees a global minimum.

We will cover three types of GD:

• Batch Gradient Descent (BGD)
• Stochastic Gradient Descent (SGD)
• Mini-batch Gradient Descent

Batch/Full Gradient Descent (BGD)

Some basic knowledge of Calculus is needed:

• Derivatives of a function with one variable using various rules (Power, Product, Quotient): https://www.youtube.com/watch?v=QqF3i1pnyzU
• Derivatives of multivariable functions: https://www.khanacademy.org/math/multivariable-calculus/multivariable-derivatives
• Non-Differentiable functions https://en.wikipedia.org/wiki/Differentiable_function

An example (see Derivative on Wikipedia for more):

• the black curve is the graph of the function $f(x) = x^3 - 3x$

• the red line is the tangent line to the function at $x = -0.6$

• the slope of the tangent line is the derivative of the function at that point: $f'(x) = 3x^2 - 3$, so $f'(-0.6) = -1.92$. The slope is negative, so the function goes down as $x$ increases — gradient descent would step to the right

  

BGD for Linear Regression

Two terms are used throughout this lecture:

• Step: one update of $\mathbf{w}$ (take one step downhill).
• Epoch: one pass over all $n$ training rows.

Parts 1-3 below make one step; part 4 repeats them.

1. Calculate partial derivative for ONE weight $w_j$

   Recall the cost function for linear regression is MSE:

   $$MSE(\mathbf{w}, b) = \frac{1}{n} \sum_{i=1}^n (\mathbf{w}^\top\mathbf{x}^{(i)} + b - y^{(i)})^2$$

   The term in the parentheses is the error of row $i$: the model's prediction $\hat{y}^{(i)} = \mathbf{w}^\top\mathbf{x}^{(i)} + b$ minus the actual value $y^{(i)}$, that is, predicted − actual. A positive error means the model predicts too high, a negative one too low. MSE squares every row's error and averages them.

   The partial derivative of the cost function with regard to ONE model parameter $w_j$ is:

   $$\frac{\partial}{\partial w_j}MSE(\mathbf{w}, b)= \frac{2}{n} \sum_{i=1}^n (\mathbf{w}^\top \mathbf{x}^{(i)} + b - y^{(i)}) x_j^{(i)}$$

   In words: multiply each row's error by that row's value of feature $j$, $x_j^{(i)}$, add up these products over all $n$ rows, and multiply the sum by $\frac{2}{n}$. The partial derivative for the bias $b$ is the same with $x_j^{(i)}$ replaced by 1, so it is just $\frac{2}{n}$ times the sum of the errors.

   This means how much the cost function is going to change if you just change $w_j$ a little bit - "what's the slope of the mountain under my feet if I face $w_j$ direction".

   NOTE: here you need to use ALL n rows to calculate the partial derivative given $w_j$ affect every row of the dataset.

   A housing example. Take the three features from the Linear Regression lecture: $x_1$ = population, $x_2$ = median_income, $x_3$ = housing_median_age. To keep the arithmetic short, use only the first two rows ($n = 2$):

   $i$	$x_1^{(i)}$ (population)	$x_2^{(i)}$ (median income)	$x_3^{(i)}$ (age)	$y^{(i)}$
   1	322	8.3252	41	452600
   2	2401	8.3014	21	358500

   Suppose gradient descent is currently at $w_1 = 0$, $w_2 = 40000$, $w_3 = 2000$, $b = 10000$. First compute each row's error, predicted − actual, where the prediction is $\hat{y}^{(i)} = \mathbf{w}^\top \mathbf{x}^{(i)} + b$:

   • row 1: predicted $= 0 \cdot 322 + 40000 \cdot 8.3252 + 2000 \cdot 41 + 10000 = 425008$, so error $= 425008 - 452600 = -27592$ (predicts too low)
   • row 2: predicted $= 0 \cdot 2401 + 40000 \cdot 8.3014 + 2000 \cdot 21 + 10000 = 384056$, so error $= 384056 - 358500 = 25556$ (predicts too high)

   Now apply the formula above for $w_2$ (median income), so $j = 2$, in three calculations:

   1. For each row, multiply the row's error by the row's median income $x_2^{(i)}$: row 1 gives $(-27592)(8.3252)$, row 2 gives $(25556)(8.3014)$.
   2. Add up these products.
   3. Multiply the sum by $\frac{2}{n}$, which is $\frac{2}{2} = 1$ here.
   $$\frac{\partial}{\partial w_2}MSE = \frac{2}{2}\big[(-27592)(8.3252) + (25556)(8.3014)\big] = -229709 + 212151 = -17558$$

2. Calculate partial derivatives for ALL parameters to get Gradient Vector $\nabla_{\mathbf{w}}$ (Nubla)

   Part 1 gives the slope for one weight. The gradient vector collects the slopes for all the parameters: every weight and the bias.

   To write all of them at once, we use a trick: append a column of 1s to $\mathbf{X}$, so the bias $b$ is absorbed as the first entry of $\mathbf{w}$ (the code below does this with X_with_bias). Then the model is just $\mathbf{X}\mathbf{w}$, and the cost is $MSE(\mathbf{w})$.

   ALL partial derivatives for the cost function can be calculated as the Gradient Vector of the cost function $\nabla_{\mathbf{w}}$ (Nubla):

   $$\nabla_{\mathbf{w}}MSE(\mathbf{w})=\begin{bmatrix} \frac{\partial}{\partial b}MSE(\mathbf{w})\\ \frac{\partial}{\partial w_1}MSE(\mathbf{w})\\ \vdots \\ \frac{\partial}{\partial w_d}MSE(\mathbf{w}) \end{bmatrix}= \frac{2}{n}\mathbf{X}^\top(\mathbf{X}\mathbf{w} - \mathbf{y})$$

   IMPORTANT NOTE: $\nabla_{\mathbf{w}}MSE(\mathbf{w})$ is calculated in EACH STEP of the Gradient Descent using the FULL training dataset.

   Therefore, Batch GD (maybe Full GD is a better name) could be very slow on very large dataset!!!

   Back to the housing example. Repeat the same three calculations from part 1 for $w_1$, $w_3$, and $b$. Only the number each error is multiplied by changes: the row's population $x_1^{(i)}$ for $w_1$, the row's age $x_3^{(i)}$ for $w_3$, and $1$ for $b$:

   $$\frac{\partial}{\partial w_1}MSE = \frac{2}{2}\big[(-27592)(322) + (25556)(2401)\big] = 52475332$$ $$\frac{\partial}{\partial w_3}MSE = \frac{2}{2}\big[(-27592)(41) + (25556)(21)\big] = -594596$$ $$\frac{\partial}{\partial b}MSE = \frac{2}{2}\big[(-27592)(1) + (25556)(1)\big] = -2036$$

   Together, the four slopes form the gradient vector of this example (the bias first, then the weights in order):

   $$\nabla MSE = \begin{bmatrix} \frac{\partial}{\partial b}MSE \\ \frac{\partial}{\partial w_1}MSE \\ \frac{\partial}{\partial w_2}MSE \\ \frac{\partial}{\partial w_3}MSE \end{bmatrix} = \begin{bmatrix} -2036 \\ 52475332 \\ -17558 \\ -594596 \end{bmatrix}$$

   The two rows share the same errors, but each weight's slope is different because each weight multiplies a different feature. The slope for $w_2$ is negative, so gradient descent will increase $w_2$. The slope for $w_1$ is huge because population values are in the thousands: row 2's large population (2401) dominates the sum. This is why features on very different scales make gradient descent struggle, and why we scale features later on this page.

   The same gradient can be built one row at a time. Each row has its own gradient, $2 \times \text{error} \times$ that row's features (the bias's "feature" is 1), which is the gradient the MSE would have if that row were the only training data. The BGD gradient is the average of the row gradients:

   $$\frac{1}{2}\left(\underbrace{2(-27592)\begin{bmatrix} 1 \\ 322 \\ 8.3252 \\ 41 \end{bmatrix}}_{\text{row 1}} + \underbrace{2(25556)\begin{bmatrix} 1 \\ 2401 \\ 8.3014 \\ 21 \end{bmatrix}}_{\text{row 2}}\right) = \frac{1}{2}\left(\begin{bmatrix} -55184 \\ -17769248 \\ -459418 \\ -2262544 \end{bmatrix} + \begin{bmatrix} 51112 \\ 122719912 \\ 424301 \\ 1073352 \end{bmatrix}\right) = \begin{bmatrix} -2036 \\ 52475332 \\ -17558 \\ -594596 \end{bmatrix}$$

   The two rows disagree: row 1 predicts too low and its gradient is all negative (increase every parameter), while row 2 predicts too high and its gradient is all positive (decrease every parameter). Averaging keeps whichever pull is stronger: row 1 wins for $b$, $w_2$, and $w_3$, row 2 wins for $w_1$. With $n$ rows, the gradient is the average of $n$ row gradients.

3. Update all parameters for the next step with learning rate Eta

   Now, we can calculate the GD step:

   $$\mathbf{w} \leftarrow \mathbf{w} - \eta \nabla_{\mathbf{w}}MSE(\mathbf{w})$$

   Here eta $\eta$ is the learning rate, which decides the size of the step. A common default value for the learning rate eta is 0.1 or 0.01. eta is often tuned via grid search from 0.1 to 10^-6.

   Updating all the weights at once is taking the steepest direction. The steepest direction just isn't along any single weight: it is a mix of all of them, with each parameter moving by $-\eta$ times its own slope.

4. Repeat parts 1-3 (one step each time) with updated $\mathbf{w}$ and full training dataset, until a stopping condition is met:

   • the maximum number of epochs is reached, or
   • the gradient is close to zero: every slope in the gradient vector, like the four in the housing example, is almost zero, so the ground is almost flat in every direction and you are at or very near the bottom. To check this with one number, compare the gradient's length, its norm $\|\nabla MSE\| = \sqrt{\text{sum of the squared slopes}}$, with a tiny tolerance $\epsilon$. In the housing example the norm is about 52,478,704, so gradient descent is far from done.

   In practice you use both: a large epoch limit as a safety net, and the tolerance to stop early once there is little left to gain. See Stopping via Tolerance below.

   Each step of batch GD uses the full training set, so in batch GD one step is one epoch, and an epoch limit is also a step limit. This stops being true for SGD, where one epoch is many steps.

   You will also see the word iteration, and it is ambiguous. Some code uses it for a step. Scikit-learn's SGD estimators use it for an epoch: in SGDRegressor, max_iter is the maximum number of epochs, n_iter_ is the number of epochs actually run, and t_ counts the steps. This lecture says step or epoch instead.

Next, we will compare linear regression and BGD and see they will generate the same optimal result:

# let's generate some linear looking data
# the ideal model apparently is y = 5 + 3x

import numpy as np
import matplotlib.pyplot as plt
plt.style.use('seaborn-v0_8')

np.random.seed(42) # fix the random seed so that each run generates the same set of random numbers
X = 2 * np.random.rand(100, 1)  # generate 100 random numbers between 0 and 2 with shape (100, 1)
y = 5 + 3 * X + np.random.randn(100, 1)  # generate 100 random numbers from a normal distribution

fig, ax = plt.subplots()
ax.plot(X, y, ".")
ax.plot(X, 5+3*X)

# Linear Regression via Scikit-Learn
import numpy as np

np.random.seed(42) # fix the random seed so that each run generates the same set of random numbers
X = 2 * np.random.rand(100, 1)  # generate 100 random numbers between 0 and 2 with shape (100, 1)
y = 5 + 3 * X + np.random.randn(100, 1)  # generate 100 random numbers from a normal distribution

from sklearn.linear_model import LinearRegression
lin_reg = LinearRegression()
lin_reg.fit(X, y) 
print(lin_reg.intercept_, lin_reg.coef_)

Output:
[5.21509616] [[2.77011339]]

try BGD:

# Batch Gradient Descent from Scratch
import numpy as np

np.random.seed(42) # fix the random seed so that each run generates the same set of random numbers
X = 2 * np.random.rand(100, 1)  # generate 100 random numbers between 0 and 2 with shape (100, 1)
y = 5 + 3 * X + np.random.randn(100, 1)  # generate 100 random numbers from a normal distribution

X_with_bias = np.c_[np.ones((100,)), X] # add bias term coefficient x0 = 1 for each instance

eta = 0.1  # the learning rate
n_epochs = 1000  # how many steps we take: in BGD, one step is one epoch
n = len(X_with_bias)  # number of training data

w = np.random.randn(2, 1) # random initialization of w with shape (2, 1): w[0] is the bias b, w[1] the weight
print(f'Initial w is \n{w}')

for epoch in range(n_epochs):  # go over the whole training set 1000 times
    gradients = 2/n * X_with_bias.T.dot(X_with_bias.dot(w) - y)  # this is gradient vector nubla
    w = w - eta * gradients  # update w

# after 1000 steps we get exactly the same result of normal equation
print(f'w after {n_epochs} epochs is \n{w}')

Output:

Initial w is 
[[0.01300189]
 [1.45353408]]
w after 1000 epochs is 
[[5.21509616]
 [2.77011339]]

Effects of Learning Rate and Number of Epochs

Next, we show how different values of learning rates and different numbers of epochs can affect BGD. The value of learning rate eta $\eta$ can be fine tuned via grid search.

• BGD converged to the optimal model

• take a smaller step, same 1000 steps - not converged to the optimal yet

• take large steps, same 1000 steps - not even converging!!

• small step and stop too early to get to the minimum

Stopping via Tolerance $\epsilon$ (epsilon)

How to set the number of epochs to avoid the situation like stopped too early?

One solution is:

1. Set a very large number of epochs

2. Set a tiny threshold called tolerance $\epsilon$ (epsilon).

3. Stop when the norm of the gradients is smaller than the tolerance.

If $\mathbf{u} \in \mathbb{R}^n$, the Norm of $\mathbf{u}$ is denoted as $\Vert \mathbf{u} \Vert$, which is defined as the length or magnitude of the vector and is calculated using $\Vert \mathbf{u} \Vert = \sqrt{u_1^2+u_2^2+\ldots + u_n^2}$.

Tolerance is often specified using scientific notation: tol=1e-3

Scikit-learn's SGDRegressor applies its tol to the loss rather than the gradient: it stops once the loss has failed to beat its best value so far by at least tol for n_iter_no_change epochs in a row (5 by default). The idea is the same: stop when another step barely helps.

# you should know scientific notation https://www.mathsisfun.com/numbers/scientific-notation.html
# 1.e-04 is 0.0001
np.format_float_scientific(0.0001)

Stochastic Gradient Descent (SGD)

BGD vs. SGD

1. Batch Gradient Descent (BGD) uses the FULL training set to calculate the gradients at each step. SGD instead only randomly pick ONE instance to compute the gradients in each step.

Think about BGD is like “look at all directions and choose sort of the average of all descents from all directions to take the next step” vs. SGD just randomly choose one direction and take the step.

2. BGD gently decreases to the optimal if the learning rate and number of epochs are appropriate, while SGD tends to bounce up and down and decreases on average.

3. BGD can become very slow when the training set is large while SGD is much faster

NOTE: an epoch is one pass over the full training set. In BGD one step is one epoch, so the number of steps and the number of epochs are the same. In SGD one epoch is $n$ steps, each on one randomly picked instance: with 100 training instances, 50 epochs means 5,000 parameter updates.

For example, with 100 training rows and 10 epochs:

Both look at the same 1,000 rows, so per epoch they do about the same amount of work. The difference is how that work is spent: BGD takes 10 careful steps, each using the average gradient of all 100 rows, while SGD takes 1,000 rough steps, each using one row's gradient. This is what "BGD is slow" means: it must go through the whole training set before taking a single step. With 10 million rows, BGD takes one step per pass over the data, while SGD has taken 10 million steps by the end of its first epoch and is often already close to a good answer.

The housing example with SGD. Go back to the two housing rows from the BGD walkthrough above and the same starting point, $b = 10000$, $w_1 = 0$, $w_2 = 40000$, $w_3 = 2000$, with learning rate $\eta = 10^{-7}$ (tiny, because the features are not scaled). With these 2 rows, one BGD epoch is 1 step: both rows go into one averaged gradient and $\mathbf{w}$ is updated once. One SGD epoch is 2 steps, each using only one row. With one row, $n = 1$, so the gradient formula becomes $2 \times \text{error} \times x_j$.

1. Step 1: pick row 1 (SGD picks at random; say row 1 comes first).

   The error is $425008 - 452600 = -27592$, so the model predicts too low. The gradient from row 1 alone, in the order $(b, w_1, w_2, w_3)$, is the row 1 gradient from the BGD walkthrough:

   $$2(-27592)\begin{bmatrix} 1 \\ 322 \\ 8.3252 \\ 41 \end{bmatrix} = \begin{bmatrix} -55184 \\ -17769248 \\ -459418 \\ -2262544 \end{bmatrix}$$

   All slopes are negative, so the update $\mathbf{w} \leftarrow \mathbf{w} - \eta \nabla$ moves every parameter up: $b = 10000.0055$, $w_1 = 1.777$, $w_2 = 40000.046$, $w_3 = 2000.226$.

2. Step 2: pick row 2, using the new $\mathbf{w}$.

   Row 2's error has grown, from 25556 to about 29828 (in BGD, row 2's error stays 25556, because every row is evaluated at the same starting $\mathbf{w}$ before any update). Raising $w_1$ helped row 1, but row 2 has a large population, so $1.777 \times 2401 \approx 4266$ was added to its prediction, which was already too high. The gradient from row 2 alone is

   $$2(29828)\begin{bmatrix} 1 \\ 2401 \\ 8.3014 \\ 21 \end{bmatrix} = \begin{bmatrix} 59655 \\ 143231822 \\ 495221 \\ 1252757 \end{bmatrix}$$

   Now every parameter moves down, and $w_1$ swings past zero: $b = 9999.9996$, $w_1 = -12.546$, $w_2 = 39999.996$, $w_3 = 2000.101$.

After one epoch, compared with one epoch of BGD from the same start and learning rate:

Three things to notice:

• Each SGD step listens to one row only. Row 1 pulls every parameter up, then row 2 pulls them back down. This tug-of-war is the bouncing you see in the SGD plots below.
• BGD averages the rows; SGD follows them one at a time. BGD steps once in the average direction of the row gradients from part 2. SGD steps in row 1's direction, then row 2's. And SGD's row 2 gradient is not the row 2 gradient from part 2: BGD evaluates every row at the same starting $\mathbf{w}$ (row 2's error 25556), while SGD reaches row 2 after step 1 has already changed $\mathbf{w}$ (error 29828).
• Feature scaling matters here too. $w_1$ (population) swings by more than 14 in one step, while $b$ and $w_2$ barely move: one row's large population dominates every step.

In summary, the final parameter values of SGD may not be optimal but often good enough given the huge gain on training speed.

# Stachastic Gradient Descent with fixed learning rate
# SGD is able to get close to the optimal
import numpy as np

n_epochs = 50 # same meaning as n_epochs in BGD: passes over the training set
eta = 0.1  # the learning rate

n = len(X_with_bias) # number of training data, 100 in this case
w = np.random.randn(2, 1) # random initialization: w[0] is the bias b, w[1] the weight

for epoch in range(n_epochs): # go over the whole training set 50 times
    for i in range(n): # note: range(stop) begins at 0 and ends at stop – 1
        random_index = np.random.randint(n)  # get a random number between 0 and n
        xi = X_with_bias[random_index:random_index+1]  # randomly choose ONE instance x_i from the training dataset
        yi = y[random_index:random_index+1]  # choose the corresponding y_i
        gradients = 2 * xi.T.dot(xi.dot(w) - yi)  # calculate the gradient for one instance
        w = w - eta * gradients  # update the gradient

print(w)

Annealing (Learning Schedule)

When getting to the minimum, in order to keep SGD from jumping too much, one solution is to gradually reduce the value of the learning rate, which is called annealing (slowly cooling down).

The intuition is to take smaller and smaller steps when getting close to the optimal.

Learning schedule refers to the function that controls the learning rate at each step.

• SGD with fixed learning rate (eta=0.1) - lots of bouncing around

  

• SGD with a learning schedule - compared with the example above, the steps became smaller and smaller and more concentrated:

  

Mini-batch Gradient Descent

SGD uses only ONE random instance to calculate gradients at each step, while BGD uses the FULL training set. Mini-BGD uses a random subset of the training data to calculate gradients at each step - less erratic than SGD but much faster than BGD.

SGD via Sklearn

Next, we use the housing dataset to train a SGD regressor and compare with the Linear regressor (RMSE: 81324) we trained in the previous lecture. If you don't have housing.csv from that lecture, download housing.csv.zip and unzip it into the folder where you run the code.

Example 1: this model's RMSE is a huge number - horrible result - the model was not even converging!

Lessons learned:

• simple regressor without tuning often lead to horrible results
• feature scaling is very important for SGD

# SGD without feature scaling
# stop after 1000 epochs, or once the loss fails to beat its best value
# by at least 0.001/1e-3 for 5 epochs in a row (n_iter_no_change=5 by default)
# eta0 is the initial eta with default invscaling learning schedule

import pandas as pd

df = pd.read_csv("housing.csv")

# for simplicity let's only choose population, median_income, housing_median_age to train the model 
# split train and test sets
X_housing = df[['population', 'median_income','housing_median_age']]
y_housing = df['median_house_value']

# split train/test
from sklearn.model_selection import train_test_split
X_train_housing, X_test_housing, y_train_housing, y_test_housing = train_test_split(X_housing, y_housing, test_size=0.2, random_state=42)

# simple sgd regressor
# this model's RMSE is a huge number - horrible result - the model was not even converging
# feature scaling and parameter tuning needed!!
from sklearn.linear_model import SGDRegressor
sgd_reg_housing = SGDRegressor(max_iter=1000, tol=1e-3, eta0=0.1)

sgd_reg_housing.fit(X_train_housing, y_train_housing)  # train the model

y_pred_housing = sgd_reg_housing.predict(X_test_housing)  # prediction

from sklearn.metrics import mean_squared_error

sgd_mse_housing = mean_squared_error(y_test_housing, y_pred_housing) 
sgd_rmse_housing = np.sqrt(sgd_mse_housing)
print(sgd_rmse_housing)

Example 2: SGD with feature scaling and parameter tuning via grid search - much better results, matching the normal equation!

import numpy as np
import pandas as pd

df = pd.read_csv("housing.csv")

# for simplicity let's only choose population, median_income, housing_median_age to train the model 
# split train and test sets
X_housing = df[['population', 'median_income','housing_median_age']]
y_housing = df['median_house_value']

# split train/test
from sklearn.model_selection import train_test_split
X_train_housing, X_test_housing, y_train_housing, y_test_housing = train_test_split(X_housing, y_housing, test_size=0.2, random_state=42)

# sgd regression pipeline with feature scaling and grid search
# RMSE is about 81323 - much better than without tuning
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
from sklearn.linear_model import SGDRegressor

sgd_reg_pipeline = Pipeline(
    steps=[
        ('scaler', StandardScaler()),
        ('sgd_reg', SGDRegressor()),
    ]
)

from sklearn.model_selection import GridSearchCV

param_grid = [
    {
        'sgd_reg__max_iter':[100000, 1000000],  # if number is too small, you will get a warning
        'sgd_reg__tol':[1e-10, 1e-3],
        'sgd_reg__eta0':[0.001, 0.01]
    }
]

# sgd regression grid search
sgd_grid_search = GridSearchCV(sgd_reg_pipeline, param_grid, cv=10)

sgd_grid_search.fit(X_train_housing, y_train_housing)

# sgd regression - best model
sgd_best = sgd_grid_search.best_estimator_
print('The best model is: \n', sgd_best)

# sgd regression - make prediction
y_pred_best_housing = sgd_best.predict(X_test_housing)

# calculate RMSE
from sklearn.metrics import mean_squared_error

sgd_rmse = np.sqrt(mean_squared_error(y_test_housing, y_pred_best_housing))
print(f'rmse: {sgd_rmse}')

# compare with the normal equation: the optimum SGD is trying to reach
from sklearn.linear_model import LinearRegression

lin_reg_housing = LinearRegression()
lin_reg_housing.fit(X_train_housing, y_train_housing)
y_pred_lin_housing = lin_reg_housing.predict(X_test_housing)
lin_rmse = np.sqrt(mean_squared_error(y_test_housing, y_pred_lin_housing))
print(f'normal equation rmse: {lin_rmse}')

The tuned SGD reaches an RMSE of about 81,323, essentially the same as the normal equation (LinearRegression), and the coefficients match too. Tuning did not make SGD better than linear regression: it made SGD converge to the same optimum the normal equation computes directly. Grid search can only get SGD as close to that optimum as possible. To do better than about 81,323, you need a different model or better features.

Model vs. Training Algorithm

In this lecture, we trained the same linear regression model three ways: the normal equation, batch GD, and SGD. All three found nearly the same parameters. This is the difference between a model and an algorithm:

• a model is what makes predictions: a function with parameters to fill in, such as $\hat{y} = \mathbf{w}^\top \mathbf{x} + b = w_1 x_1 + \dots + w_d x_d + b$
• a training algorithm is how the parameters get filled in from data, such as the normal equation or gradient descent
• a loss function defines what "best" parameters means, such as MSE

Models and training algorithms mix and match:

In scikit-learn, each class (an estimator) bundles these choices. LinearRegression is a linear model with MSE loss, trained by solving the least-squares problem directly. The SGDRegressor we used above is named after its algorithm, but it is also a linear model, trained by SGD. Its loss and penalty parameters choose the rest. With the defaults (loss='squared_error', penalty='l2'), it is linear regression with a small amount of regularization added. SGDRegressor(penalty=None) is plain linear regression trained by SGD.

Appendix: Plotting Code

The concept figures at the top of the page:

import numpy as np
import matplotlib.pyplot as plt
plt.style.use('seaborn-v0_8-white')


def gradient_descent(grad, w, eta, n_steps):
    """Run gradient descent on one parameter and return every w it visits."""
    path = [w]
    for _ in range(n_steps):
        w = w - eta * grad(w)
        path.append(w)
    return np.array(path)


def draw_steps(ax, cost, path, color='tab:blue', curve=-0.3):
    """Draw each step as a dot on the cost curve with an arrow to the next one."""
    ax.scatter(path, cost(path), s=90, color=color, edgecolor='black', zorder=3)
    for a, b in zip(path[:-1], path[1:]):
        ax.annotate('', xy=(b, cost(b)), xytext=(a, cost(a)),
                    arrowprops=dict(arrowstyle='->', color='gray', connectionstyle=f'arc3,rad={curve}'))


def style_axes(ax):
    """A concept figure: label the axes but hide the tick values."""
    ax.set_xticks([])
    ax.set_yticks([])
    ax.set_xlabel(r'$w$', fontsize=14)
    ax.set_ylabel('Cost', fontsize=12)


# a bowl-shaped (convex) cost function with its minimum at w = 3
def cost(w):
    return (w - 3) ** 2 + 1


def cost_grad(w):
    return 2 * (w - 3)


w_plot = np.linspace(-0.2, 6.2, 300)

# figure 1: gradient descent steps toward the minimum
path = gradient_descent(cost_grad, w=0.2, eta=0.15, n_steps=8)
fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(w_plot, cost(w_plot), color='black', linewidth=2)
draw_steps(ax, cost, path)
ax.scatter([3], [cost(3)], s=140, color='gold', edgecolor='black', zorder=4)
ax.annotate('random initial value', xy=(path[0], cost(path[0])), xytext=(0.6, 9.5),
            arrowprops=dict(arrowstyle='->', color='black'))
ax.annotate('learning step', xy=((path[1] + path[2]) / 2, cost(path[1]) + 0.3), xytext=(2.0, 6.5),
            arrowprops=dict(arrowstyle='->', color='black'))
ax.annotate('minimum', xy=(3, cost(3)), xytext=(3.6, 3.2), arrowprops=dict(arrowstyle='->', color='black'))
style_axes(ax)

# figure 2: a learning rate that is too small - many tiny steps, still far from the minimum
path = gradient_descent(cost_grad, w=0.2, eta=0.02, n_steps=12)
fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(w_plot, cost(w_plot), color='black', linewidth=2)
draw_steps(ax, cost, path)
ax.set_title(r'learning rate too small ($\eta$ = 0.02): 12 steps and still far from the minimum')
style_axes(ax)

# figure 3: a learning rate that is too large - each step overshoots and the cost grows
path = gradient_descent(cost_grad, w=2.4, eta=1.1, n_steps=6)
fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(w_plot, cost(w_plot), color='black', linewidth=2)
draw_steps(ax, cost, path, color='tab:red', curve=0)
ax.set_title(r'learning rate too large ($\eta$ = 1.1): each step overshoots and the cost grows')
style_axes(ax)


# figure 4: a cost function that is not convex
def bumpy_cost(w):
    return (1.7 - 0.9 * np.exp(-(w - 1.5) ** 2 / 0.5)       # local minimum
            - 1.6 * np.exp(-(w - 4.5) ** 2 / 0.7)           # global minimum
            + 0.9 / (1 + np.exp(-(w - 6.2) * 2.0))          # step up to the plateau
            + 0.6 * np.exp(-w * 1.2)                        # steep wall on the left
            + 0.25 * np.maximum(w - 9, 0) ** 2)             # rising right edge


def bumpy_grad(w, h=1e-5):
    return (bumpy_cost(w + h) - bumpy_cost(w - h)) / (2 * h)


w_plot = np.linspace(0.2, 10.6, 800)
fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(w_plot, bumpy_cost(w_plot), color='black', linewidth=2)
# start on the left: gradient descent gets stuck in the local minimum
draw_steps(ax, bumpy_cost, gradient_descent(bumpy_grad, w=0.35, eta=0.3, n_steps=6))
# start on the right: the slope flattens out, so the steps get tiny on the plateau
draw_steps(ax, bumpy_cost, gradient_descent(bumpy_grad, w=10.5, eta=1.0, n_steps=6), curve=0.4)
for w, label in [(1.53, 'local minimum'), (4.49, 'global minimum')]:
    ax.vlines(w, -0.05, bumpy_cost(w), linestyle='--', color='gray')
    ax.text(w, -0.3, label, ha='center', va='bottom', fontsize=11)
ax.text(8.2, bumpy_cost(8.2) + 0.2, 'plateau', ha='center', fontsize=11)
ax.set_ylim(-0.35, 3.5)
style_axes(ax)


# figure 5: the derivative is the slope of the tangent line
def f(x):
    return x ** 3 - 3 * x


def f_prime(x):
    return 3 * x ** 2 - 3


x_plot = np.linspace(-2.2, 2.2, 300)
x0 = -0.6
fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(x_plot, f(x_plot), color='black', linewidth=2, label=r'function $f(x) = x^3 - 3x$')
ax.plot(x_plot, f(x0) + f_prime(x0) * (x_plot - x0), color='tab:red', linewidth=2,
        label=f'tangent line at x = {x0}: slope = f\'({x0}) = {f_prime(x0):.2f}')
ax.scatter([x0], [f(x0)], s=90, color='tab:red', edgecolor='black', zorder=3)
ax.set_ylim(-4, 4)
ax.set_xlabel('x')
ax.set_ylabel('f(x)')
ax.legend(loc='upper center')

The learning-path figures:

# Plotting code for the learning-path figures above.
# Self-contained: it regenerates the same data (seed 42) as the page code.
import numpy as np
import matplotlib.pyplot as plt
plt.style.use('seaborn-v0_8')

def make_data():
    np.random.seed(42)  # same data and same initial w as the BGD example
    X = 2 * np.random.rand(100, 1)
    y = 5 + 3 * X + np.random.randn(100, 1)
    X_with_bias = np.c_[np.ones((100,)), X]
    w = np.random.randn(2, 1)  # initial w [[0.01300189], [1.45353408]]: w[0] is the bias b, w[1] the weight
    return X, y, X_with_bias, w

def plot_path(ws, X, y, title):
    # left: (b, w1) after each step; right: the data and the lines from the first 10 steps
    ws = np.array(ws)
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 4))
    ax1.plot(ws[:, 0], ws[:, 1], '.')
    ax1.set_xlabel(r'$b$')
    ax1.set_ylabel(r'$w_1$')
    ax1.set_title(title, fontsize=15)
    ax2.plot(X, y, 'd', color='tan', markeredgecolor='gray', markersize=4)
    for b, w1 in ws[:10]:
        ax2.plot(X, b + w1 * X)
    ax2.set_title('first 10 steps', fontsize=15)
    return fig

def bgd(eta, n_epochs):
    X, y, X_with_bias, w = make_data()
    n = len(X_with_bias)
    ws = [w.ravel()]
    for epoch in range(n_epochs):  # in BGD, one step per epoch
        gradients = 2/n * X_with_bias.T.dot(X_with_bias.dot(w) - y)
        w = w - eta * gradients
        ws.append(w.ravel())
    return plot_path(ws, X, y, f'eta={eta}, epochs={n_epochs}')

def sgd(n_epochs, schedule):
    X, y, X_with_bias, w = make_data()
    n = len(X_with_bias)
    ws = [w.ravel()]
    for epoch in range(n_epochs):
        for i in range(n):
            eta = schedule(epoch * n + i)
            random_index = np.random.randint(n)
            xi = X_with_bias[random_index:random_index+1]
            yi = y[random_index:random_index+1]
            gradients = 2 * xi.T.dot(xi.dot(w) - yi)
            w = w - eta * gradients
            ws.append(w.ravel())
    return ws, X, y, eta

# Batch GD with different learning rates and numbers of epochs
bgd(eta=0.1, n_epochs=1000)   # converged to the optimal model
bgd(eta=0.01, n_epochs=1000)  # smaller steps: not converged yet
bgd(eta=0.5, n_epochs=1000)   # steps too large: diverges
bgd(eta=0.01, n_epochs=100)   # small steps and stopped too early

# SGD with a fixed learning rate
ws, X, y, eta = sgd(n_epochs=50, schedule=lambda t: 0.1)
plot_path(ws, X, y, 'eta=0.1, epochs=50')

# SGD with a learning schedule: eta = t0 / (t + t1) shrinks at every step
t0, t1 = 5, 50
ws, X, y, eta = sgd(n_epochs=50, schedule=lambda t: t0 / (t + t1))
plot_path(ws, X, y, f'learning schedule (final eta={eta:.5f}), epochs=50')