Ordinary least squares derivation

Overview

Ordinary Least Squares (OLS) is the analytical method for finding the optimal parameters in linear regression by minimizing the sum of squared residuals. Unlike gradient descent (iterative), OLS gives us a closed-form solution through calculus.

Think: "We punish the model quadratically for being wrong, so it tries hardest to fix big mistakes."

The Setup: What Are We Optimizing?

Problem Statement

Given:

• Training data: $\{(x_1, y_1), (x_2, y_2), \ldots, (x_m, y_m)\}$ where $x_i \in \mathbb{R}^n$, $y_i \in \mathbb{R}$
• Model: $\hat{y}_i = \theta^T x_i = \theta_0 + \theta_1 x_{i1} + \cdots + \theta_n x_{in}$

Find: Parameter vector $\theta \in \mathbb{R}^{n+1}$ that minimizes the cost function:

$J(\theta) = \frac{1}{2m} \sum_{i=1}^{m} (h_\theta(x_i) - y_i)^2$

Why the $\frac{1}{2m}$?

• The $m$ averages over examples (makes cost comparable across dataset sizes)
• The $\frac{1}{2}$ cancels with the derivative of $x^2$ (convenience, doesn't change the optimal $\theta$)

Matrix Formulation

Stack everything into matrices (this is WHERE the power of linear algebra enters):

1 & x_{11} & x_{12} & \cdots & x_{1n} \\ 1 & x_{21} & x_{22} & \cdots & x_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ 1 & x_{m1} & x_{m2} & \cdots & x_{mn} \end{bmatrix}_{m \times (n+1)}, \quad y = \begin{bmatrix} y_1 \\ y_2 \\ \vdots \ y_m \end{bmatrix}_{m \times 1}, \quad \theta = \begin{bmatrix} \theta_0 \\ \theta_1 \\ \vdots \\ \theta_n \end{bmatrix}_{(n+1) \times 1}$$ **Why the column of1s?** That's the bias term $\theta_0$ (intercept). By adding it to $X$, we handle bias and weights uniformly: $\hat{y} = X\theta$. Now our cost becomes: $$J(\theta) = \frac{1}{2m} \|X\theta - y\|^2 = \frac{1}{2m}(X\theta - y)^T(X\theta - y)$$ ## Derivation from First Principles ### Step 1: Expand the Squared Error $$J(\theta) = \frac{1}{2m}(X\theta - y)^T(X\theta - y)$$ Let's expand this using $(A - B)^T(A - B) = A^TA - A^TB - B^TA + B^TB$: $$J(\theta) = \frac{1}{2m}\left[(X\theta)^T(X\theta) - (X\theta)^Ty - y^T(X\theta) + y^Ty\right]$$ **Why this step?** We need to see all the $\theta$ terms explicitly to differentiate. Since $(X\theta)^Ty$ is a scalar, it equals its transpose $y^T(X\theta)$: $$J(\theta) = \frac{1}{2m}\left[\theta^TX\theta - 2y^TX\theta + y^Ty\right]$$ ### Step 2: Take the Gradient To minimize, we need $\nabla_\theta J(\theta) = 0$. > [!definition] Matrix Calculus Rules We Need > 1. $\frac{\partial}{\partial \theta}(a^T\theta) = a$ (linear term) > 2. $\frac{\partial}{\partial \theta}(\theta^TA\theta) = 2A\theta$ (quadratic term, when $A$ is symmetric) > 3. $\frac{\partial}{\partial \theta}(c) = 0$ (constant term) Applying these: $$\nabla_\theta J(\theta) = \frac{1}{2m}\left[2X^TX\theta - 2X^Ty + 0\right] = \frac{1}{m}\left[X^TX\theta - X^Ty\right]$$ **Why does the2 appear?** The quadratic form $\theta^TX\theta$ differentiates to $2X^TX\theta$ (both sides contribute symmetrically). ### Step 3: Set Gradient to Zero $$\frac{1}{m}\left[X^TX\theta - X^Ty\right] = 0$$ $$X^TX\theta = X^Ty$$ > [!formula] The Normal Equation > $$\boxed{\theta^* = (X^TX)^{-1}X^Ty}$$ This is the **analytical solution** to linear regression. The matrix $(X^TX)^{-1}X^T$ is called the ==Moore-Penrose pseudoinverse==. **WHY "Normal"?** In geometry, "normal" means perpendicular. At the optimum, the residual vector $(X\theta^* - y)$ is **orthogonal** to the column space of $X$. That's why $X^T(X\theta^* - y) = 0$. ## Worked Examples > [!example] Example 1: Simple Linear Regression (1D) > **Setup**: Fit $y = \theta_0 + \theta_1 x$ to data: $(1, 2), (2, 4), (3, 5)$ **Step 1**: Build matrices $$X = \begin{bmatrix} 1 & 1 \\ 1 & 2 \\ 1 & 3 \end{bmatrix}, \quad y = \begin{bmatrix} 2 \\ 4 \\ 5 \end{bmatrix}$$ **Step 2**: Compute $X^TX$ $$X^TX = \begin{bmatrix} 1 & 1 & 1 \\ 1 & 2 & 3 \end{bmatrix} \begin{bmatrix} 1 & 1 \\ 1 & 2 \\ 1 & 3 \end{bmatrix} = \begin{bmatrix} 3 & 6 \\ 6 & 14 \end{bmatrix}$$ **Why this step?** We need to invert $X^TX$ for the normal equation. **Step 3**: Compute $X^Ty$ $$X^Ty = \begin{bmatrix} 1 & 1 \\ 1 & 2 & 3 \end{bmatrix} \begin{bmatrix} 2 \\ 4 \\ 5 \end{bmatrix} = \begin{bmatrix} 11 \ 25 \end{bmatrix}$$ **Step 4**: Invert $X^TX$ $$\text{det}(X^TX) = 3(14) - 6(6) = 42 - 36 = 6$$ $$(X^TX)^{-1} = \frac{1}{6}\begin{bmatrix} 14 & -6 \\ -6 & 3 \end{bmatrix} = \begin{bmatrix} 7/3 & -1 \\ -1 & 1/2 \end{bmatrix}$$ **Step 5**: Solve $$\theta^* = \begin{bmatrix} 7/3 & -1 \ -1 & 1/2 \end{bmatrix} \begin{bmatrix} 11 \\ 25 \end{bmatrix} = \begin{bmatrix} 77/3 - 25 \\ -11 + 25/2 \end{bmatrix} = \begin{bmatrix} 2/3 \\ 3/2 \end{bmatrix}$$ > **Answer**: $\theta_0 = 2/3$, $\theta_1 = 3/2$, so $\hat{y} = \frac{2}{3} + \frac{3}{2}x$ > > **Verification**: Plug in $x=2$: $\hat{y} = 2/3 + 3 = 3.67$ (close to actual $y=4$ ✓) > > [!example] Example 2: Why (X^TX) Might Not Be Invertible > **Scenario**: 3 data points, 2 features, but feature 2 = $2 \times$ feature 1 (==multicollinearity==) > $$X = \begin{bmatrix} 1 & 1 & 2 \\ 1 & 2 & 4 \\ 1 & 3 & 6 \end{bmatrix}$$ **Compute**: $$X^TX = \begin{bmatrix} 3 & 6 & 12 \\ 6 & 14 & 28 \\ 12 & 28 & 56 \end{bmatrix}$$ **Check determinant**: Row 3 =2 × Row 2, so the matrix is **singular** (det = 0). **Why this matters**: When features are linearly dependent, infinitely many $\theta$ minimize the cost. The problem is ==ill-posed==. **Solution**: Use ==regularization== (Ridge/Lasso) or remove redundant features. ## Common Mistakes > [!mistake] Mistake 1: Forgetting the Bias Column > **Wrong thinking**: "I'll just use $X$ with my features, the algorithm will figure out the intercept." **Why it feels right**: We write $y = mx + b$ separately in basic algebra. **The fix**: Linear algebra formulation requires the bias as an explicit column of1s. Without it, your line is **forced through the origin** (underfitting). **Test**: Fit $y = 2x + 3$ data without bias. You'll get $\theta_1 \approx 2.something$ but predictions will be off by ~3. > [!mistake] Mistake 2: Using OLS When $m < n$ (Underdetermined System) > **Wrong thinking**: "More features = better model, let me add100 features for10 data points." **Why it feels right**: More parameters should capture more patterns. **Reality**: When $m < n+1$, $X^TX$ is singular (not enough equations to uniquely solve for all $\theta$). You get ==overfitting== — the model memorizes noise. **The fix**: Use regularization or reduce features to $n < m$. > [!mistake] Mistake 3: Assuming OLS is Always Faster Than Gradient Descent > **Wrong thinking**: "Closed form = instant solution." **Reality**: Inverting $X^TX$ is $O(n^3)$ and storing it is $O(n^2)$. For $n > 10,000$ features, this is **slower and more memory-intensive** than iterative methods. **Rule of thumb**: Use OLS when $n< 1000$ and data fits in memory. Otherwise, use gradient descent or stochastic variants. ## When the Normal Equation Fails $(X^TX)$ is invertible **if and only if** $X$ has ==full column rank== (all columns are linearly independent). **Failure cases**: 1. **Redundant features**: Two features are perfectly corelated (e.g., temperature in °C and °F) 2. **More features than samples**: $n > m$ (underdetermined) 3. **Duplicate rows**: Same data point repeated **Detection**: Compute `np.linalg.matrix_rank(X.T @ X)`. If rank $< n+1$, you have a problem. **Fixes**: - Remove corelated features (check correlation matrix) - Add ==regularization==: $\theta^* = (X^TX + \lambda I)^{-1}X^Ty$ (Ridge regression) - Use dimensionality reduction (PCA) ## Computational Considerations ### Complexity Analysis | Operation | Complexity | Why? | |-----------|--------|---| | $X^TX$ | $O(mn^2)$ | Multiply $m \times n$ with $n \times m$ | | $(X^TX)^{-1}$ | $O(n^3)$ | Gaussian elimination or LU decomposition | | $(X^TX)^{-1}X^Ty$ | $O(n^2m)$ | Matrix-vector multiplications | | **Total** | $O(mn^2 + n^3)$ | Dominated by $n^3$ when $n$ is large | **Comparison with Gradient Descent**: - GD: $O(kmn)$ where $k$ = iterations (typically $k \ll n$) - GD wins when $n > 1000$ or data doesn't fit in RAM ### Numerical Stability > [!intuition] Why Direct Inversion Can Fail > Even when $(X^TX)$ is theoretically invertible, floating-point arithmetic can make it **numerically singular** if: > - Features have vastly different scales (e.g., age in years vs. income in dollars) > - Condition number of $X^TX$ is high (matrix is "nearly" singular) **Better approach**: Use ==QR decomposition== or ==SVD== instead of direct inversion. Libraries like `numpy.linalg.lstsq` do this automatically. ```python # Bad: theta = np.linalg.inv(X.T @ X) @ X.T @ y # Good: theta = np.linalg.lstsq(X, y, rcond=None)[0] ``` ## Geometric Interpretation The solution $\theta^*$ makes the prediction vector $X\theta^*$ the **orthogonal projection** of $y$ onto the column space of $X$. **Why?** The residual $e = y - X\theta^*$ is perpendicular to every column of $X$: $$X^T e = X^T(y - X\theta^*) = 0$$ This is exactly our normal equation! Geometrically, we're finding the **closest point** in the column space of $X$ to the target $y$. > [!mnemonic] "PLOT" for OLS Steps > - **P**roblem: Set up $X$, $y$ (don't forget bias!) > - **L**inear algebra: Compute $X^TX$ and $X^Ty$ > - **O**ptimize: Solve $(X^TX)\theta = X^Ty$ > - **T**est: Check invertibility and validate on data ## Connections - [[Linear Regression Fundamentals]] — the model this solves - [[Gradient Descent]] — iterative alternative to OLS - [[Ridge Regression]] — adds $\lambda I$ to fix non-invertibility - [[Maximum Likelihood Estimation]] — OLS is MLE under Gaussian noise - [[QR Decomposition]] — numerically stable way to solve normal equations - [[Pseudoinverse]] — generalizes matrix inverse for singular matrices - [[Feature Scaling]] — prevents numerical issues in $X^TX$ - [[Bias-Variance Tradeoff]] — OLS has zero bias but can have high variance --- > [!recall]- Explain to a 12-Year-Old > Imagine you're trying to draw the "best-fit" straight line through some dots on a graph. How do you know when your line is the "best"? One way is to measure how far each dot is from your line (these are called "errors" or "mistakes"). If a dot is above the line, you measure up; if below, you measure down. But here's the problem: some mistakes are positive (above) and some are negative (below), so they might cancel out even if your line is terrible! So instead, we **square** all the mistakes. Now every mistake is positive (because negative × negative = positive). We add up all these squared mistakes and try to make that total as small as possible. That's the "least squares" part! The "ordinary" part just means we're using normal math (not anything fancy). And here's the cool part: there's a magic formula that gives you the perfect line **instantly** without guessing. It's like having the answer key before the test! We use some matrix math (fancy multiplication with tables of numbers) to find the exact slope and intercept that make the total squared error the smallest possible. That formula is $(X^TX)^{-1}X^Ty$ — it looks scary, but it's just: "Mix your data in a special way and out pops the perfect answer!" The reason this works is that we're using calculus to find where the error is at its minimum — like finding the bottom of a valley. At that bottom, if you take one tiny step in any direction, the error starts going up. That's how we know we've found the best line. --- #flashcards/ai-ml What is the Ordinary Least Squares (OLS) method? :: An analytical approach to find optimal linear regression parameters by minimizing the sum of squared residuals, yielding a closed-form solution $\theta^* = (X^TX)^{-1}X^Ty$ through calculus. Why do we square the residuals in OLS instead of using absolute values? ::: (1) Squaring eliminates sign cancellation, (2) penalizes large errors more heavily (outlier sensitivity), (3) results in differentiable cost function for calculus, and (4) corresponds to maximum likelihood estimation under Gaussian noise. What is the Normal Equation? ::: $\theta^* = (X^TX)^{-1}X^Ty$, the closed-form solution to linear regression obtained by setting the gradient of the cost function to zero. Called "normal" because the residual vector is orthogonal (normal) to the column space of X. Why do we add a column of 1s to the design matrix X? ::: To incorporate the bias term (intercept) $\theta_0$ into the matrix formulation, allowing uniform treatment of bias and weights in the equation $\hat{y} = X\theta$. When is $(X^TX)$ not invertible? ::: When X does not have full column rank: (1) redundant/correlated features, (2) more features than samples ($n > m$), or (3) duplicate data points. This makes the system underdetermined with infinitely many solutions. What is the computational complexity of solving OLS with the Normal Equation? ::: $O(mn^2 + n^3)$, dominated by the $O(n^3)$ matrix inversion step. Computing $X^TX$ takes $O(mn^2)$ and the final multiplication takes $O(n^2m)$. When should you use gradient descent instead of OLS? ::: When $n > 1000$ features (inversion becomes expensive), data doesn't fit in memory, or you need online learning. Gradient descent has complexity $O(kmn)$ where $k$ is iterations, typically much faster for large $n$. What does it mean geometrically that $X^T(y - X\theta^*) = 0$? ::: The residual vector $e = y - X\theta^*$ is orthogonal to the column space of X. The prediction $X\theta^*$ is the orthogonal projection of y onto the column space of X — the closest point in that space to y. Why might directly computing $(X^TX)^{-1}$ fail numerically even when it's theoretically invertible? ::: Features at different scales or high condition number can make the matrix numerically singular due to floating-point precision limits. Better to use QR decomposition or SVD (e.g., `np.linalg.lstsq`). How does Ridge regression fix the non-invertibility problem? ::: Adds regularization term: $\theta^* = (X^TX + \lambda I)^{-1}X^Ty$. The $\lambda I$ term ensures $(X^TX + \lambda I)$ is always invertible (positive definite) even when $X^TX$ is singular. What is multicollinearity and why does it break OLS? ::: When features are linearly dependent (one can be written as a combination of others), making $X^TX$ singular. This creates infinitely many optimal solutions, making the problem ill-posed. Derive why $\frac{\partial}{\partial \theta}(\theta^TX\theta) = 2X^TX\theta$ ::: Let $A = X^TX$ (symmetric). Expand: $\theta^TA\theta = \sum_i \sum_j \theta_i A_{ij} \theta_j$. Taking $\frac{\partial}{\partial \theta_k}$ gives $\sum_j A_{kj}\theta_j + \sum_i A_{ik}\theta_i = 2\sum_j A_{kj}\theta_j = 2(A\theta)_k$ since $A$ is symmetric. Thus $\nabla_\theta(\theta^TA\theta) = 2A\theta$. ## 🖼️ Concept Map ```mermaid flowchart TD OLS[Ordinary Least Squares] -->|minimizes| COST[Cost function J theta] COST -->|sum of| SQ[Squared residuals] SQ -->|avoids sign cancel| WHY[Squaring rationale] WHY -->|equals| MLE[Max likelihood under Gaussian] COST -->|written as| MAT[Matrix form X theta minus y] MAT -->|adds| BIAS[Column of 1s for intercept] MAT -->|expand terms| EXP[Expanded quadratic] EXP -->|set gradient zero| GRAD[Gradient of J] GRAD -->|solves to| NEQ[Normal equations] NEQ -->|closed-form| SOL[theta equals XtX inv Xt y] OLS -->|alternative to| GD[Gradient descent iterative] SOL -->|optimal params for| LR[Linear regression] ``` ## 🔊 Hinglish (regional understanding) > [!intuition]- Hinglish mein samjho > Linear regression mein sabse important chez hai **best fit line** dhoondhna. Lekin ye line "best" kaise decide karein? OLS (Ordinary Least Squares) ka matlab hai ki hum sabhi **errors ko square karke add** karte hain, aur jo line is total ko sabse **chhota** bana de, wahi best hai. > > **Kyon square karein?** Agar simple differencelete (y - ŷ) toh positive aur negative errors cancel ho jate. Square karne se sab positive ho jaate hain, aur bade errors ko **zyada punishment** milti hai (outliers handle karne ke liye). Plus, ye **calculus mein smooth** rehta hai, toh differentiate karke minimum dhoondhna easy ho jata hai. > > **Normal Equation** ye magic formula hai: θ* = (X^TX)^(-1)X^Ty. Ye ek **closed-form solution** hai matlab ek baar mein seedha answer mil jata hai, gradient descent ki tarah baar-baar iterate nahi karna padta. Par ek problem: agar **features zyada hain** (n > 1000) ya **data bahut bada** hai, toh ye inverse nikalna **slow** ho jata hai (O(n³) complexity). Tab gradient descent better hai. > > Practical tip: Hamesha X mein **bias column** (column of 1s) add karo taki intercept (θ₀) bhi solve ho jaaye. Aur agar (X^TX) ** ![[audio/2.2.03-Ordinary-least-squares-derivation.mp3]]