Concepts 🚀 Project 🧮 Labs

Linear Regression maths — the equation behind .fit()

model.fit() is not magic. It solves y = Xβ + ε using least squares. Understand the equation, the cost function, and the geometry.

script.py v3.1

import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression

# ── The equation: y = β₀ + β₁x₁ + β₂x₂ + ... + ε ──
# β = coefficients (what we learn)
# ε = error (what we can't explain)

# ── Generate data with KNOWN coefficients ──
np.random.seed(42)
n = 100
x = np.random.uniform(0, 10, n)
# TRUE equation: y = 2 + 3x + noise
true_intercept = 2
true_slope = 3
noise = np.random.normal(0, 2, n)  # ε
y = true_intercept + true_slope * x + noise

print('TRUE equation: y = 2 + 3x + noise')
print(f'TRUE coefficients: intercept={true_intercept}, slope={true_slope}')

# ── Fit LinearRegression — it should RECOVER the true coefficients ──
model = LinearRegression()
model.fit(x.reshape(-1, 1), y)
print(f'\nLEARNED coefficients: intercept={model.intercept_:.2f}, slope={model.coef_[0]:.2f}')
print(f'Close to true? intercept diff={abs(model.intercept_ - true_intercept):.2f}, slope diff={abs(model.coef_[0] - true_slope):.2f}')

# ════════════════════════════════════════════════════════════════════
#  THE MATHS: how does .fit() find the coefficients?
# ════════════════════════════════════════════════════════════════════
print('\n' + '='*60)
print('THE MATHS: Least Squares')
print('='*60)
print('Goal: find β that minimizes Σ(y_actual - y_predicted)²')
print('This is the "cost function" or "loss function".')
print()
print('Matrix form: y = Xβ + ε')
print('  y = column vector of targets (n×1)')
print('  X = matrix of features (n×p) — first column is all 1s for the intercept')
print('  β = column vector of coefficients (p×1)')
print('  ε = error vector (n×1)')
print()
print('Solution (the NORMAL EQUATION):')
print('  β = (XᵀX)⁻¹Xᵀy')
print('  This is what .fit() computes — in one matrix operation.')

# ── Verify by computing β manually ──
X_matrix = np.column_stack([np.ones(n), x])  # add intercept column
beta_manual = np.linalg.inv(X_matrix.T @ X_matrix) @ X_matrix.T @ y
print(f'\nManual computation: β = {beta_manual}')
print(f'sklearn:            intercept={model.intercept_:.4f}, slope={model.coef_[0]:.4f}')
print('SAME! sklearn just uses a more numerically stable algorithm (SVD).')

# ════════════════════════════════════════════════════════════════════
#  THE COST FUNCTION — what we're minimizing
# ════════════════════════════════════════════════════════════════════
print('\n' + '='*60)
print('THE COST FUNCTION: Mean Squared Error')
print('='*60)
print('MSE(β) = (1/n) * Σ(yᵢ - β₀ - β₁xᵢ)²')
print('This is a BOWL shape (convex) — has one minimum.')
print('The minimum is where the gradient = 0.')
print('Gradient descent: take steps in the -gradient direction.')

# Visualize the cost function for different slopes
slopes = np.linspace(0, 6, 100)
mse_values = [np.mean((y - true_intercept - s*x)**2) for s in slopes]
fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(slopes, mse_values, color='#a0c000', linewidth=2)
ax.axvline(true_slope, color='green', linestyle='--', label=f'True slope={true_slope}')
ax.axvline(model.coef_[0], color='red', linestyle='--', label=f'Learned slope={model.coef_[0]:.2f}')
ax.set_xlabel('Slope (β₁)')
ax.set_ylabel('MSE (cost)')
ax.set_title('Cost Function — the bowl shape')
ax.legend()
plt.tight_layout()
plt.savefig('cost_function.png', dpi=100, bbox_inches='tight')
plt.show()
print('\nThe BOWL shape: minimum at the true slope (3). sklearn found it (2.97).')

# ════════════════════════════════════════════════════════════════════
#  GEOMETRY — projection onto a column space
# ════════════════════════════════════════════════════════════════════
print('\n' + '='*60)
print('GEOMETRY: Linear Regression as Projection')
print('='*60)
print('Think of X as defining a PLANE (column space).')
print('y is a point NOT on the plane (because of noise ε).')
print('The PREDICTION ŷ = Xβ is the PROJECTION of y onto the plane.')
print('The residual (y - ŷ) is perpendicular to the plane.')
print('This is why it\'s called "least squares" — minimizes the squared perpendicular distance.')

Python output click Run to execute

Click "Run" to execute the code on the left.
The output will appear here.

🧠 How it works

Linear Regression maths: EQUATION: y = β₀ + β₁x₁ + ... + ε. β are coefficients, ε is error. COST FUNCTION: MSE = (1/n)Σ(y - ŷ)². Convex bowl shape — one minimum. SOLUTION (Normal Equation): β = (XᵀX)⁻¹Xᵀy. Closed-form — computed in one matrix op. XᵀX is the Gram matrix. Xᵀy is the feature-target correlation. Inverse makes it a 'divide'. sklearn uses SVD (more stable than inverse). GEOMETRY: X defines a plane. y is a point above the plane. ŷ is the projection. Residual (y - ŷ) is perpendicular to the plane. Least squares = minimize squared perpendicular distance. GRADIENT DESCENT (alternative): start at random β, take steps in -gradient direction. Used when XᵀX is too big to invert (millions of features). ASSUMPTIONS (often violated, often don't matter): Linearity, independence, homoscedasticity (constant variance), normal residuals. KEY INSIGHT: .fit() is not magic. It's matrix algebra. Understand it = debug models.

✏️ Try these changes

Increase noise from 2 to 10
→ learned coefficients will be further from true — more error
Change true_slope from 3 to -2
→ negative relationship — cost function still works
Add a non-linear term: y = 2 + 3x + 0.5x² + noise
→ LinearRegression will underfit — need PolynomialFeatures
Compute R² manually: 1 - SS_res/SS_tot
→ matches model.score()

← Previous lab

EDA visualization — 6 plots that tell you everything

Next lab →

Gradient Descent — the engine behind neural networks