Objective: Derive and code the gradient descent update rule for linear regression by hand, train it on synthetic data, and confirm it converges to the same coefficients scikit-learn finds analytically. All in runnable, in-browser Python.
Setup
Every code block below runs live in your browser via Pyodide. Click Run to execute real
Python and see real output. We'll generate a synthetic single-feature regression dataset with
sklearn.datasets.make_regression, so no external data or installation is required.
Task 1: Generate Synthetic Data
Control the noise level yourself so it's easy to sanity-check whether your implementation is working later.
import numpy as np
from sklearn.datasets import make_regression
X, y = make_regression(n_samples=200, n_features=1, noise=15, random_state=42)
X = X.flatten() # single feature, work with a 1D array for clarity
print("X shape:", X.shape)
print("y shape:", y.shape)
print("First 5 X values:", np.round(X[:5], 3))
print("First 5 y values:", np.round(y[:5], 3))
Task 2: Implement the MSE Cost Function
The model we're fitting is the simple line y_pred = w * X + b. The cost is mean squared error,
averaged across all n examples.
import numpy as np
from sklearn.datasets import make_regression
X, y = make_regression(n_samples=200, n_features=1, noise=15, random_state=42)
X = X.flatten()
def compute_cost(X, y, w, b):
y_pred = w * X + b
return np.mean((y_pred - y) ** 2)
# Sanity check: cost at w=0, b=0 should just be the variance-like average of y^2
print("Cost at w=0, b=0:", round(compute_cost(X, y, 0.0, 0.0), 2))
print("Cost at w=1, b=0:", round(compute_cost(X, y, 1.0, 0.0), 2))
Task 3: Implement Gradient Descent by Hand
Taking the partial derivative of the MSE cost with respect to w and b gives the
update directions below. If you get stuck on the math: the cost is a function of y_pred, and
y_pred is a function of w and b.
Differentiating
(w*X + b - y)^2 with respect to w brings down a factor of X via the
chain rule; differentiating with respect to b does not. That's the only structural difference
between dw and db below. Run this block to train the model and print the loss every
100 iterations.
import numpy as np
from sklearn.datasets import make_regression
X, y = make_regression(n_samples=200, n_features=1, noise=15, random_state=42)
X = X.flatten()
def compute_cost(X, y, w, b):
y_pred = w * X + b
return np.mean((y_pred - y) ** 2)
def compute_gradients(X, y, w, b):
n = len(X)
y_pred = w * X + b
dw = (2 / n) * np.sum((y_pred - y) * X)
db = (2 / n) * np.sum(y_pred - y)
return dw, db
def train(X, y, lr=0.1, n_iterations=1000):
w, b = 0.0, 0.0
cost_history = []
for i in range(n_iterations):
cost = compute_cost(X, y, w, b)
cost_history.append(cost)
if i % 100 == 0:
print(f"Iteration {i:4d} | cost = {cost:.3f} | w = {w:.4f} | b = {b:.4f}")
dw, db = compute_gradients(X, y, w, b)
w -= lr * dw
b -= lr * db
return w, b, cost_history
w, b, cost_history = train(X, y, lr=0.1, n_iterations=1000)
print()
print(f"Final learned w = {w:.4f}, b = {b:.4f}")
A common stumbling block is the learning rate: too high and the cost diverges or oscillates instead of
decreasing; too low and it barely moves in 1000 iterations. If your printed cost explodes to huge numbers,
halve lr and re-run. Since make_regression output isn't tiny in scale, 0.01-0.1 is
usually a reasonable starting range.
Task 4: Plot the Loss Curve
Plotting cost against iteration should show a steep drop in the first handful of steps, followed by a long, flat tail as the parameters converge.
import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import make_regression
X, y = make_regression(n_samples=200, n_features=1, noise=15, random_state=42)
X = X.flatten()
def compute_cost(X, y, w, b):
y_pred = w * X + b
return np.mean((y_pred - y) ** 2)
def compute_gradients(X, y, w, b):
n = len(X)
y_pred = w * X + b
dw = (2 / n) * np.sum((y_pred - y) * X)
db = (2 / n) * np.sum(y_pred - y)
return dw, db
def train(X, y, lr=0.1, n_iterations=1000):
w, b = 0.0, 0.0
cost_history = []
for i in range(n_iterations):
cost_history.append(compute_cost(X, y, w, b))
dw, db = compute_gradients(X, y, w, b)
w -= lr * dw
b -= lr * db
return w, b, cost_history
w, b, cost_history = train(X, y, lr=0.1, n_iterations=1000)
plt.plot(cost_history)
plt.xlabel("Iteration")
plt.ylabel("MSE")
plt.title("Gradient descent loss curve")
show_plot()
Task 5: Compare Against scikit-learn
The moment of truth: fit scikit-learn's analytic solver on the exact same X and y and
print both sets of parameters side by side.
import numpy as np
from sklearn.datasets import make_regression
from sklearn.linear_model import LinearRegression
X, y = make_regression(n_samples=200, n_features=1, noise=15, random_state=42)
X = X.flatten()
def compute_cost(X, y, w, b):
y_pred = w * X + b
return np.mean((y_pred - y) ** 2)
def compute_gradients(X, y, w, b):
n = len(X)
y_pred = w * X + b
dw = (2 / n) * np.sum((y_pred - y) * X)
db = (2 / n) * np.sum(y_pred - y)
return dw, db
def train(X, y, lr=0.1, n_iterations=1000):
w, b = 0.0, 0.0
for i in range(n_iterations):
dw, db = compute_gradients(X, y, w, b)
w -= lr * dw
b -= lr * db
return w, b
w_gd, b_gd = train(X, y, lr=0.1, n_iterations=1000)
reg = LinearRegression()
reg.fit(X.reshape(-1, 1), y)
w_sklearn, b_sklearn = reg.coef_[0], reg.intercept_
print(f"{'':12s} {'slope (w)':>12s} {'intercept (b)':>15s}")
print(f"{'Gradient descent':12s} {w_gd:12.4f} {b_gd:15.4f}")
print(f"{'scikit-learn':12s} {w_sklearn:12.4f} {b_sklearn:15.4f}")
print()
print(f"Difference in slope: {abs(w_gd - w_sklearn):.5f}")
print(f"Difference in intercept: {abs(b_gd - b_sklearn):.5f}")
If your gradient descent converged, w and b should land within a small fraction of a
percent of scikit-learn's values. They will almost never be bit-for-bit identical, since gradient descent is an
iterative approximation while scikit-learn solves the normal equations directly. If your numbers are way off,
the most likely culprits are an unstable learning rate, too few iterations, or a sign error in one of the two
gradient formulas. Double-check the chain rule derivation in Task 3.