3. Simple Linear Regression

3. Simple Linear Regression — From Intuition to OLS and Gradient Descent

Detailed notes based on 2 simpleLinearRegression.ipynb. They preserve the notebook's example, derive the formulas it uses, explain every major code section, and replace fragile or misleading practices with safer alternatives.

1. Learning goals

You should finish able to:

  • interpret the slope, intercept, prediction, and residual;
  • derive the closed-form OLS slope and intercept;
  • explain why the fitted line passes through $(\bar x,\bar y)$;
  • write the design-matrix form and explain why direct inversion is fragile;
  • derive the gradient of MSE with respect to $m$ and $b$;
  • perform one gradient-descent update by hand;
  • explain the learning rate and feature scaling;
  • distinguish Batch GD, SGD, and mini-batch GD;
  • reconcile the notebook's full-data OLS and training-data GD values;
  • implement both approaches with fully commented Python.

2. The model: what it says

Simple linear regression uses one feature $x$ to model the conditional mean of a numerical target $Y$:

$$ E[Y\mid X=x]=\beta_0+\beta_1x. $$

For an individual observation:

$$ y_i=\beta_0+\beta_1x_i+\varepsilon_i, $$

and the fitted prediction is:

$$ \hat y_i=b+mx_i. $$

Symbol Name Meaning
$x_i$ feature observed input for row $i$
$y_i$ target observed output
$\hat y_i$ prediction point on the fitted line
$m$ or $\hat\beta_1$ slope estimated target change per one-unit $x$ increase
$b$ or $\hat\beta_0$ intercept predicted target at $x=0$
$\varepsilon_i$ population error unobserved deviation from the population line
$e_i$ residual observed deviation $y_i-\hat y_i$

What does the slope mean?

$$ m=\frac{\Delta\hat y}{\Delta x}. $$

If $m=2.92$, increasing $x$ by one unit changes the prediction by $2.92$ units:

$$ \hat y(x+1)-\hat y(x)=m. $$

The units of $m$ are "target units per feature unit." If $x$ is floor area and $y$ is price, the slope is price per area unit.

What does the intercept mean?

$$ b=\hat y(0). $$

It is the line's vertical-axis crossing. It has a practical interpretation only if $x=0$ is meaningful and not unreasonable extrapolation.

Why "simple"?

"Simple" means one predictor, not easy data, a small dataset, or an unimportant model.

3. Geometric intuition

Every candidate pair $(m,b)$ defines a line. The best-fit rule compares those lines using the vertical gaps between each observed point and its fitted point.

flowchart TD A["Choose candidate m and b"] --> B["Compute ŷᵢ = mxᵢ + b"] B --> C["Compute residuals eᵢ = yᵢ - ŷᵢ"] C --> D["Square and aggregate residuals"] D --> E{"Is this the minimum?"} E -->|"No"| F["Change m and b"] F --> B E -->|"Yes"| G["Best-fitting line"] classDef parameter fill:#ede9fe,stroke:#7c3aed,color:#3b0764,stroke-width:2px classDef calculate fill:#dbeafe,stroke:#2563eb,color:#1e3a8a classDef cost fill:#fef3c7,stroke:#d97706,color:#78350f,stroke-width:2px classDef update fill:#fee2e2,stroke:#dc2626,color:#7f1d1d classDef finish fill:#dcfce7,stroke:#16a34a,color:#14532d,stroke-width:2px class A parameter class B,C calculate class D,E cost class F update class G finish

For one observation:

$$ e_i=y_i-(mx_i+b). $$

A point above the line has $e_i>0$: the prediction was too small. A point below the line has $e_i<0$: the prediction was too large.

Why vertical, not shortest, distances?

Ordinary regression treats $x$ as given and models randomness in $y$. If both axes contain substantial measurement error, the usual slope can be biased and orthogonal-distance or errors-in-variables methods may be preferable.

Parameter-space intuition

The data plot lives in $(x,y)$-space. Optimization lives in $(m,b,J)$-space. Each point on the loss surface represents an entire candidate line.

flowchart LR A["Data space: points and lines"] --> B["One line ↔ one pair (m, b)"] B --> C["Parameter space: loss J(m, b)"] C --> D["Lowest point ↔ best line"] classDef data fill:#cffafe,stroke:#0891b2,color:#164e63,stroke-width:2px classDef map fill:#fef3c7,stroke:#d97706,color:#78350f,stroke-width:2px classDef loss fill:#fce7f3,stroke:#db2777,color:#831843,stroke-width:2px classDef best fill:#dcfce7,stroke:#16a34a,color:#14532d,stroke-width:2px class A data class B map class C loss class D best

4. Residuals and the cost function

Sum of squared errors

$$ \operatorname{SSE}(m,b)=\sum_{i=1}^{n}[y_i-(mx_i+b)]^2. $$

Mean squared error

$$ \operatorname{MSE}(m,b)=\frac1n\sum_{i=1}^{n}[y_i-(mx_i+b)]^2. $$

Half mean squared error

$$ J(m,b)=\frac{1}{2n}\sum_{i=1}^{n}[y_i-(mx_i+b)]^2. $$

All three have the same minimizer because they differ only by positive constants. The factor $1/2$ is a calculus convenience: differentiating the square produces a $2$, which then cancels.

This distinction matters in code. The notebook's custom class calculates MSE with $1/n$, so its gradient correctly includes a factor of $2$. The markdown uses $1/(2n)$, whose gradient omits that factor. Both conventions are valid when loss and gradient are consistent.

Why squared loss?

  • Positive and negative residuals cannot cancel.
  • A residual of 10 contributes $100$, while a residual of 1 contributes $1$.
  • The objective is differentiable.
  • The linear-regression loss is convex.

The sensitivity to large errors is both a feature and a weakness. If extreme observations are common or data errors occur, compare robust alternatives such as Huber or absolute-error regression.

5. Deriving ordinary least squares

We minimize:

$$ J(m,b)=\sum_{i=1}^{n}[y_i-(mx_i+b)]^2. $$

The positive scaling constant is omitted because it cannot change the minimizer.

Step 1: differentiate with respect to $b$

$$ \frac{\partial J}{\partial b}=-2\sum_i[y_i-(mx_i+b)]. $$

Set it to zero:

$$ \sum_i y_i-m\sum_i x_i-nb=0. $$

Divide by $n$:

$$ \bar y-m\bar x-b=0, $$

so:

$$ \boxed{b=\bar y-m\bar x}. $$

This immediately proves that the fitted OLS line passes through $(\bar x,\bar y)$:

$$ m\bar x+b=m\bar x+\bar y-m\bar x=\bar y. $$

Step 2: differentiate with respect to $m$

$$ \frac{\partial J}{\partial m}=-2\sum_i x_i[y_i-(mx_i+b)]. $$

After substituting $b=\bar y-m\bar x$ and collecting centered terms:

$$ \boxed{m=\frac{\sum_i(x_i-\bar x)(y_i-\bar y)}{\sum_i(x_i-\bar x)^2}}. $$

Then:

$$ \boxed{b=\bar y-m\bar x}. $$

Statistical interpretation of the slope

Using sample covariance and variance with matching denominators:

$$ m=\frac{\operatorname{Cov}(X,Y)}{\operatorname{Var}(X)}=r_{xy}\frac{s_y}{s_x}. $$

Consequences:

  • If covariance is positive, slope is positive.
  • If correlation is zero, the simple OLS slope is zero.
  • Rescaling $x$ changes the numeric slope.
  • If all $x_i$ are equal, the denominator is zero and slope cannot be identified.

Worked example

Let:

$$ x=[1,2,3],\qquad y=[2,3,5]. $$

Then:

$$ \bar x=2,\qquad \bar y=\frac{10}{3}. $$

$$ \sum(x_i-\bar x)(y_i-\bar y)=(-1)\left(-\frac43\right)+0+\left(\frac53\right)=3. $$

$$ \sum(x_i-\bar x)^2=1+0+1=2. $$

Therefore:

$$ m=\frac32=1.5,\qquad b=\frac{10}{3}-1.5(2)=\frac13. $$

The fitted line is:

$$ \hat y=\frac13+1.5x. $$

Predictions are approximately $[1.833,3.333,4.833]$; residuals are $[0.167,-0.333,0.167]$. Notice that residuals sum to zero.

6. Matrix OLS and numerical stability

For multiple features, include a ones column in the design matrix:

$$ \tilde X=\begin{bmatrix}1 & x_{11} & \cdots & x_{1p}\\ 1 & x_{21} & \cdots & x_{2p}\\ \vdots & \vdots & \ddots & \vdots\\ 1 & x_{n1} & \cdots & x_{np}\end{bmatrix}. $$

The model is:

$$ \hat y=\tilde X\boldsymbol{\beta}. $$

The least-squares objective is:

$$ \lVert y-\tilde X\boldsymbol{\beta}\rVert_2^2. $$

Differentiating:

$$ \nabla_{\boldsymbol{\beta}}J=-2\tilde X^\top(y-\tilde X\boldsymbol{\beta}). $$

Setting it to zero gives the normal equations:

$$ \tilde X^\top\tilde X\hat{\boldsymbol{\beta}}=\tilde X^\top y. $$

If invertible:

$$ \hat{\boldsymbol{\beta}}=(\tilde X^\top\tilde X)^{-1}\tilde X^\top y. $$

Why the notebook's direct inverse is educational but fragile

The notebook uses:

beta = np.linalg.inv(X_bias.T @ X_bias) @ X_bias.T @ y

It mirrors the textbook formula, but production code should avoid forming the inverse because:

  • floating-point roundoff is amplified;
  • $X^\top X$ squares the condition number;
  • exact multicollinearity makes it singular;
  • near-collinearity makes coefficients unstable.

Prefer:

beta, residual_sums, rank, singular_values = np.linalg.lstsq(
    X_bias,
    y,
    rcond=None,
)

Scikit-learn's dense LinearRegression wraps a stable least-squares solver rather than naïvely evaluating the explicit inverse.

Projection intuition

OLS selects $\hat y$ as the projection of $y$ onto the column space of $X$. At the optimum:

$$ X^\top e=0. $$

Thus residuals are orthogonal to every included feature column. If an intercept column is included:

$$ \mathbf{1}^\top e=\sum_i e_i=0. $$

flowchart TD A["Target vector y"] --> B["Project onto column space of X"] B --> C["Fitted vector ŷ = Xβ̂"] A --> D["Residual vector e = y - ŷ"] C --> E["Orthogonality: Xᵀe = 0"] D --> E classDef target fill:#fce7f3,stroke:#db2777,color:#831843,stroke-width:2px classDef operation fill:#fef3c7,stroke:#d97706,color:#78350f classDef fit fill:#dbeafe,stroke:#2563eb,color:#1e3a8a,stroke-width:2px classDef residual fill:#fee2e2,stroke:#dc2626,color:#7f1d1d classDef result fill:#dcfce7,stroke:#16a34a,color:#14532d,stroke-width:2px class A target class B operation class C fit class D residual class E result

7. Gradient descent from first principles

Gradient descent reaches the same OLS minimum iteratively.

Using MSE:

$$ J(m,b)=\frac1n\sum_i[y_i-(mx_i+b)]^2. $$

Define:

$$ \hat y_i=mx_i+b,\qquad e_i=y_i-\hat y_i. $$

The derivatives are:

$$ \frac{\partial J}{\partial m}=-\frac{2}{n}\sum_i x_i(y_i-\hat y_i), $$

$$ \frac{\partial J}{\partial b}=-\frac{2}{n}\sum_i(y_i-\hat y_i). $$

Updates:

$$ m_{\text{new}}=m_{\text{old}}-\alpha\frac{\partial J}{\partial m}, $$

$$ b_{\text{new}}=b_{\text{old}}-\alpha\frac{\partial J}{\partial b}. $$

Directional intuition

The gradient points in the local direction of fastest increase. Subtracting it moves downhill.

flowchart TD A["Initialize m and b"] --> B["Predict all rows"] B --> C["Measure MSE"] C --> D["Calculate dm and db"] D --> E["Update opposite the gradient"] E --> F{"Converged or iteration limit?"} F -->|"No"| B F -->|"Yes"| G["Return fitted parameters"] classDef init fill:#ede9fe,stroke:#7c3aed,color:#3b0764,stroke-width:2px classDef compute fill:#dbeafe,stroke:#2563eb,color:#1e3a8a classDef gradient fill:#fef3c7,stroke:#d97706,color:#78350f,stroke-width:2px classDef update fill:#fce7f3,stroke:#db2777,color:#831843 classDef decision fill:#fee2e2,stroke:#dc2626,color:#7f1d1d classDef done fill:#dcfce7,stroke:#16a34a,color:#14532d,stroke-width:2px class A init class B,C compute class D gradient class E update class F decision class G done

One update by hand

Use the worked data $x=[1,2,3]$, $y=[2,3,5]$, start $m=b=0$, and use notebook-style MSE.

$$ \hat y=[0,0,0],\qquad e=[2,3,5]. $$

$$ \frac{\partial J}{\partial m}=-\frac23(1\cdot2+2\cdot3+3\cdot5)=-\frac{46}{3}\approx-15.333. $$

$$ \frac{\partial J}{\partial b}=-\frac23(2+3+5)=-\frac{20}{3}\approx-6.667. $$

At $\alpha=0.01$:

$$ m_{\text{new}}=0-0.01(-15.333)=0.1533, $$

$$ b_{\text{new}}=0-0.01(-6.667)=0.0667. $$

Both move toward the optimum $m=1.5,b=1/3$.

Learning-rate behavior

flowchart TD A["Learning rate α"] --> B["Too large"] A --> C["Suitable"] A --> D["Very small"] B --> B1["Overshoot, oscillate, or diverge"] C --> C1["Stable useful convergence"] D --> D1["Stable but many iterations"] classDef root fill:#ede9fe,stroke:#7c3aed,color:#3b0764,stroke-width:2px classDef bad fill:#fee2e2,stroke:#dc2626,color:#7f1d1d,stroke-width:2px classDef good fill:#dcfce7,stroke:#16a34a,color:#14532d,stroke-width:2px classDef slow fill:#fef3c7,stroke:#d97706,color:#78350f,stroke-width:2px class A root class B,B1 bad class C,C1 good class D,D1 slow

The notebook diagram says a low learning rate "may get stuck." For convex linear least squares, "very slow" is the accurate intuition; local minima are not the issue.

Why feature scaling helps GD

If one feature ranges from $0$ to $1$ and another from $0$ to $10^6$, the loss surface becomes elongated. A step size safe in the steep direction is inefficient in the flat direction. Standardization:

$$ z_j=\frac{x_j-\mu_j}{s_j} $$

improves conditioning and allows more balanced updates. Fit $\mu_j$ and $s_j$ using training data only.

Stopping criteria

Common rules:

  • absolute loss change below a tolerance;
  • gradient norm below a tolerance;
  • parameter movement below a tolerance;
  • validation loss stops improving;
  • maximum iterations reached.

A tiny loss change can occur because the learning rate is tiny, so combining criteria is safer.

8. Batch, stochastic, and mini-batch GD

Aspect Batch GD SGD Mini-batch GD
Rows per update All $n$ One Small batch
Path Smooth Noisy Moderately noisy
Update cost High Low Moderate
Vectorization Excellent Weak Excellent
Streaming Poor Excellent Possible
Common use Smaller convex problems Online learning Large modern training

The notebook table labels Batch GD "high accuracy" and SGD "low/noisy." Noise describes the optimization path, not necessarily final predictive accuracy. With suitable schedules and enough updates, SGD can approximate the same optimum and sometimes generalize well.

An epoch means one full pass through the training data. Batch GD performs one update per epoch; SGD performs approximately $n$; mini-batch GD performs approximately $n/B$ updates for batch size $B$.

9. OLS versus gradient descent

flowchart TD A["Need linear least squares"] --> B{"Small or medium dense problem?"} B -->|"Yes"| C["Use a stable OLS solver"] B -->|"No"| D{"Streaming, sparse, or extremely large?"} D -->|"Yes"| E["Use iterative or stochastic optimization"] D -->|"No"| F["Benchmark solver choices"] C --> G["No learning-rate tuning"] E --> H["Scale features and monitor convergence"] F --> I["Compare time, memory, and validation error"] classDef root fill:#ede9fe,stroke:#7c3aed,color:#3b0764,stroke-width:2px classDef decision fill:#fef3c7,stroke:#d97706,color:#78350f,stroke-width:2px classDef ols fill:#dbeafe,stroke:#2563eb,color:#1e3a8a,stroke-width:2px classDef gd fill:#fce7f3,stroke:#db2777,color:#831843,stroke-width:2px classDef result fill:#dcfce7,stroke:#16a34a,color:#14532d class A root class B,D decision class C,G ols class E,H gd class F,I result
Question Stable OLS solver Gradient descent
Exact finite-step linear algebra? Yes, subject to floating point No, iterative approximation
Learning rate? No Yes
Feature scaling essential? Not for predictions, though conditioning matters Strongly recommended
Huge/streaming data? Often less convenient Well suited
Rank deficiency? SVD/lstsq can return a solution Convergence and solution depend on setup
Convex global optimum? Yes Yes, with suitable optimization

For the notebook's 100 rows and one feature, a stable OLS solver is the natural practical choice. The custom GD class is valuable for learning how optimization works.

10. Notebook code walkthrough

10.1 Synthetic data

np.random.seed(42)
X = np.random.rand(100, 1) * 10
y = 3 * X.squeeze() + 5 + np.random.randn(100) * 2
  • The seed makes the exact random sequence reproducible.
  • $X$ has 100 rows and one feature.
  • squeeze() changes $(100,1)$ to $(100,)$.
  • The hidden slope is 3 and intercept is 5.
  • Noise prevents a perfect line and makes estimation realistic.

10.2 Scalar OLS function

The notebook's ols_simple follows the centered formulas exactly. A valuable guard would be:

if np.isclose(denominator, 0):
    raise ValueError("Slope is undefined because x has no variation.")

10.3 Matrix OLS function

np.c_[np.ones(...), X] adds an intercept column. The order of returned coefficients is $[\text{intercept},\text{slope}]$. Replace np.linalg.inv with np.linalg.lstsq for numerical stability.

10.4 Train/test split

The custom GD model fits only X_train and y_train. Therefore it should be compared with OLS fitted on the same training data, not the earlier OLS result fitted on all 100 rows.

10.5 Gradient class

Important operations:

y_pred = np.dot(X, self.weights) + self.bias

This vectorizes $\hat y=Xw+b$.

dw = -(2/n_samples) * np.dot(X.T, (y - y_pred))
db = -(2/n_samples) * np.sum(y - y_pred)

These are the gradients of MSE. The updates subtract them.

One subtlety: the original class stores history before each update. After the final update, self.weights is one step ahead of the final stored history item. That difference is tiny here, but history should ideally record a consistent state.

11. Improved commented implementation

11.1 Stable OLS functions