Linear Regression

Background

Galton (1886) coined regression studying hereditary height; Legendre formalized least squares in 1805
Foundational supervised learning algorithm for continuous outputs
The closed-form solution (normal equations) makes it unique among ML algorithms

The Concept

Linear regression fits the line (or hyperplane) through data that minimizes the sum of squared vertical distances from each point to the line. It learns one weight per input feature plus a bias term.

How it Links

Foundation for logistic regression, regularized regression (ridge/lasso), and GLMs
Normal equations are a special case of the projection theorem in linear algebra
The bias-variance trade-off is most clearly illustrated here via ridge and lasso regularization

Concrete Example

Predicting a house price from square footage: the model learns that each extra square foot adds $200, plus a base price of $50,000. Residuals are the difference between predicted and actual prices - you minimize the sum of those squared.

Prerequisites

Math you need: linear algebra basics, Bias, variance, over/under-fitting, Evaluation metrics

Learning outcomes

Derive the OLS normal equations from the mean-squared-error objective.
State the five classical assumptions and the Gauss-Markov theorem.
Recognize OLS as MLE under Gaussian noise, and ridge as the MAP estimate under a Gaussian prior.
Compare ridge, lasso, and elastic net by objective, closed form, sparsity, and behaviour under multicollinearity.
Implement OLS and ridge from scratch in NumPy, with sklearn, and as a PyTorch autograd loop.

TL;DR

Linear regression fits the hyperplane $\hat{y} = X\beta$ that minimizes the squared residual norm. The closed-form OLS solution $\hat{\beta} = (X^{\top}X)^{-1} X^{\top} y$ is the orthogonal projection of $y$ onto the column space of $X$ , and Gauss-Markov says it is the best linear unbiased estimator (BLUE) when errors are uncorrelated and homoskedastic [1]. Ridge [4], lasso [3], and elastic net [5] add penalties that trade a little bias for a lot of variance reduction, with ridge giving a closed form, lasso giving sparsity, and elastic net combining both.

Intuition

A line through a cloud of points is the simplest predictive model that is still useful. Linear regression chooses the line (or hyperplane in higher dimensions) that minimizes the sum of squared vertical residuals. Two equivalent stories explain why squared error is the right loss:

Probabilistic. If we assume the noise is i.i.d. Gaussian, $y_i = x_i^{\top}\beta + \varepsilon_i$ with $\varepsilon_i \sim \mathcal{N}(0, \sigma^2)$ , then maximum-likelihood estimation of $\beta$ is exactly minimization of the squared residual norm [2].
Geometric. Squared error is the squared Euclidean distance, so OLS finds the orthogonal projection of $y$ onto the column space of $X$ . The residual vector is the component of $y$ that lies orthogonal to every feature column.

The right picture

OLS is just orthogonal projection. The fitted vector $\hat{y}$ lives in $\operatorname{col}(X)$ and the residual $r = y - \hat{y}$ is perpendicular to every column of $X$ . That single fact gives you the normal equations, Gauss-Markov, the geometry of multicollinearity, and why ridge “shrinks” along the small-singular-value directions.

OLS as orthogonal projection of y onto col(X)

col(X) - hyperplane spanned by columns of X

x_1x_2yy_hat = X beta_hatr = y - y_hat

Normal equations: X^T r = 0 ⇔ X^T (y - X beta_hat) = 0 ⇔ beta_hat = (X^T X)^(-1) X^T y

Figure 1. The OLS fit y_hat is the orthogonal projection of y onto col(X). The residual r is perpendicular to every column of X, which is exactly the normal equations X^T r = 0.

Math

Math: derivation of the normal equations and the OLS estimator

Setup

Stack $n$ observations in $X \in \mathbb{R}^{n \times d}$ and $y \in \mathbb{R}^{n}$ . We fit $\hat{y} = X\beta$ with $\beta \in \mathbb{R}^{d}$ (an intercept is absorbed by adding a column of ones to $X$ ).

Squared-error objective

Define the residual sum of squares

\mathcal{L}(\beta) = \tfrac{1}{2}\,\lVert y - X\beta \rVert_{2}^{2} = \tfrac{1}{2}\,(y - X\beta)^{\top}(y - X\beta).

Expanding,

\mathcal{L}(\beta) = \tfrac{1}{2}\,y^{\top}y - \beta^{\top} X^{\top} y + \tfrac{1}{2}\,\beta^{\top} X^{\top} X \beta.

Setting the gradient to zero

Differentiate with respect to $\beta$ :

\nabla_{\beta} \mathcal{L}(\beta) = -\,X^{\top} y + X^{\top} X \beta.

The first-order condition $\nabla_{\beta} \mathcal{L}(\beta) = 0$ gives the normal equations

X^{\top} X \, \hat{\beta} = X^{\top} y.

If $X$ has full column rank, $X^{\top} X$ is invertible and

\boxed{\,\hat{\beta} = (X^{\top} X)^{-1} X^{\top} y.\,}

The Hessian is $\nabla^{2}_{\beta} \mathcal{L}(\beta) = X^{\top} X \succeq 0$ , so $\mathcal{L}$ is convex and $\hat{\beta}$ is a global minimum.

Geometric interpretation

The fitted vector $\hat{y} = X\hat{\beta} = X (X^{\top}X)^{-1} X^{\top} y = P_{X}\, y$ , where $P_X$ is the orthogonal projection matrix onto $\operatorname{col}(X)$ . The residual

r = y - \hat{y} = (I - P_X)\, y

lies in $\operatorname{col}(X)^{\perp}$ , so $X^{\top} r = 0$ . That is exactly the normal equations: residuals are orthogonal to every feature column.

Connection to MLE

If $y_i = x_i^{\top}\beta + \varepsilon_i$ with $\varepsilon_i \sim \mathcal{N}(0, \sigma^2)$ i.i.d., the log-likelihood is

\log p(y \mid X, \beta, \sigma^2) = -\frac{n}{2}\log(2\pi\sigma^2) - \frac{1}{2\sigma^2}\,\lVert y - X\beta \rVert_{2}^{2}.

Maximizing in $\beta$ is identical to minimizing $\lVert y - X\beta \rVert_{2}^{2}$ , so OLS is the MLE under Gaussian noise [2].

Gradient descent variant

When $X^{\top}X$ is too large to invert, run gradient descent with step size $\eta$ :

\beta_{t+1} = \beta_{t} - \eta\, \bigl(X^{\top} X \beta_{t} - X^{\top} y\bigr) = \beta_{t} - \eta\, X^{\top}(X\beta_{t} - y).

Convergence requires $\eta < 2 / \lambda_{\max}(X^{\top}X)$ . Stochastic gradient descent uses a single example or mini-batch per step.

Gauss-Markov theorem: OLS is BLUE

Statement

Assume:

Linearity: $y = X\beta + \varepsilon$ with $X$ non-random and full column rank.
Zero-mean errors: $\mathbb{E}[\varepsilon] = 0$ .
Spherical errors: $\operatorname{Cov}(\varepsilon) = \sigma^{2} I$ (uncorrelated and homoskedastic).

Then among all linear unbiased estimators of $\beta$ , the OLS estimator $\hat{\beta} = (X^{\top}X)^{-1} X^{\top} y$ has the smallest variance: it is the Best Linear Unbiased Estimator (BLUE).

Proof sketch

Any linear estimator can be written $\tilde{\beta} = C y$ for some matrix $C \in \mathbb{R}^{d \times n}$ . Write $C = (X^{\top}X)^{-1} X^{\top} + D$ , where $D$ is the deviation from the OLS row-operator.

Unbiasedness $\mathbb{E}[\tilde{\beta}] = \beta$ requires $C X = I$ , which forces $D X = 0$ .

Compute the covariance:

\operatorname{Cov}(\tilde{\beta}) = \sigma^{2} C C^{\top} = \sigma^{2}\Bigl[(X^{\top}X)^{-1} + D D^{\top}\Bigr],

using $D X = 0$ to kill the cross terms. Since $D D^{\top} \succeq 0$ ,

\operatorname{Cov}(\tilde{\beta}) - \operatorname{Cov}(\hat{\beta}_{\mathrm{OLS}}) = \sigma^{2} D D^{\top} \succeq 0,

so $\hat{\beta}_{\mathrm{OLS}}$ has the smallest variance among linear unbiased estimators. Equality requires $D = 0$ , i.e. uniqueness up to multiplicity. Note that biased estimators (ridge!) can do strictly better in MSE; Gauss-Markov is a statement about the unbiased class only [1].

The classical assumptions

OLS is consistent and well-behaved when these five assumptions hold. Real data often violates one or more, and each violation has a textbook fix.

Assumption	What it means	Common diagnostic	Standard fix
Linearity	E[y \| x] is linear in beta.	Residuals-vs-fitted plot shows no curvature.	Add polynomial features, splines, or use a non-linear model.
Independence	Observations (or errors) are independent.	Durbin-Watson on residuals; check for time-series or cluster structure.	GLS, mixed-effects, or block-resampling.
Homoskedasticity	Var(epsilon \| x) is constant.	Breusch-Pagan test; funnel-shaped residual plot.	Weighted least squares or robust (sandwich) standard errors.
Normality of residuals	epsilon is Gaussian. Needed for finite-sample inference, not consistency.	QQ-plot of residuals; Shapiro-Wilk on small n.	Bootstrap confidence intervals; do not bother for large n by CLT.
No multicollinearity	X has full column rank; columns are not near-collinear.	VIF (variance inflation factor) above 5-10; condition number of X above 30.	Drop a redundant feature, ridge regression, or PCA pre-processing.

Inference vs prediction

Strict normality matters for finite-sample t and F statistics. For prediction under a large $n$ , the central limit theorem rescues you and you can ignore non-Gaussian residuals. Always know which job your model is doing.

Regularization

Adding a penalty to the OLS objective trades a bit of bias for a lot of variance reduction, especially when features are correlated or when $d$ is close to $n$ .

Ridge regression: closed-form derivation

The ridge objective is

\mathcal{L}_{\mathrm{ridge}}(\beta) = \tfrac{1}{2}\,\lVert y - X\beta \rVert_{2}^{2} + \tfrac{\lambda}{2}\,\lVert \beta \rVert_{2}^{2}, \qquad \lambda > 0.

Setting $\nabla_{\beta}\mathcal{L}_{\mathrm{ridge}}(\beta) = -X^{\top}y + X^{\top}X\beta + \lambda \beta = 0$ ,

\boxed{\,\hat{\beta}_{\mathrm{ridge}} = (X^{\top}X + \lambda I)^{-1} X^{\top} y.\,}

The penalty term makes $X^{\top}X + \lambda I$ strictly positive definite, so the inverse always exists, even when $X^{\top}X$ is singular (e.g. $d > n$ or perfectly collinear columns). This is exactly why ridge dominates OLS under multicollinearity [4].

MAP interpretation

Place a Gaussian prior $\beta \sim \mathcal{N}(0, \tau^{2} I)$ . The posterior mode under the Gaussian likelihood is

\arg\max_{\beta}\bigl[\log p(y \mid X, \beta) + \log p(\beta)\bigr] = \arg\min_{\beta}\,\Bigl[\tfrac{1}{2\sigma^{2}}\lVert y - X\beta \rVert^{2} + \tfrac{1}{2\tau^{2}}\lVert \beta \rVert^{2}\Bigr],

with $\lambda = \sigma^{2}/\tau^{2}$ . Ridge is MAP under a Gaussian prior on the weights. Smaller $\tau^{2}$ pulls $\hat{\beta}$ harder towards zero [2].

SVD view

Let $X = U \Sigma V^{\top}$ be the thin SVD with singular values $s_1 \ge s_2 \ge \dots \ge s_d \ge 0$ . Then

\hat{\beta}_{\mathrm{ridge}} = V\,\operatorname{diag}\!\left(\frac{s_j}{s_j^{2} + \lambda}\right) U^{\top} y, \qquad \hat{\beta}_{\mathrm{OLS}} = V\,\operatorname{diag}\!\left(\frac{1}{s_j}\right) U^{\top} y.

Ridge applies the multiplicative shrinkage $s_j^{2}/(s_j^{2}+\lambda) < 1$ to the OLS coefficients, with the heaviest shrinkage on directions with small singular values. Those are exactly the directions where OLS has highest variance, so ridge buys variance reduction precisely where OLS is least stable [1].

Lasso and elastic net

The lasso replaces the squared penalty with an absolute-value penalty:

\mathcal{L}_{\mathrm{lasso}}(\beta) = \tfrac{1}{2}\,\lVert y - X\beta \rVert_{2}^{2} + \lambda \,\lVert \beta \rVert_{1}.

There is no closed form, but the objective is still convex and is solved by coordinate descent or proximal gradient. Each coordinate update reduces to the soft-thresholding operator $\operatorname{sign}(z)\,\max(|z| - \lambda, 0)$ , which sets coefficients to zero exactly when their unpenalized residual signal is below $\lambda$ . That is why lasso produces sparse solutions: it is the penalty geometry, not a heuristic [3].

The elastic net is the convex combination

\mathcal{L}_{\mathrm{enet}}(\beta) = \tfrac{1}{2}\,\lVert y - X\beta \rVert_{2}^{2} + \lambda\!\Bigl[(1-\alpha)\,\tfrac{1}{2}\lVert \beta \rVert_{2}^{2} + \alpha\,\lVert \beta \rVert_{1}\Bigr],

with $\alpha \in [0, 1]$ . It inherits lasso’s sparsity but, unlike lasso, can keep correlated groups of features together rather than arbitrarily picking one of them [5].

Method	Objective	Closed form?	Sparsity	Multicollinearity	Hyperparameter
OLS	\|\|y - X beta\|\|^2	yes (X^T X must be invertible)	no	unstable, exploding variance	none
Ridge	\|\|y - X beta\|\|^2 + lambda * \|\|beta\|\|^2	yes: (X^T X + lambda I)^(-1) X^T y	no (shrinks but does not zero)	stable, shrinks correlated features together	lambda
Lasso	\|\|y - X beta\|\|^2 + lambda * \|\|beta\|\|_1	no (coordinate descent / proximal)	yes (sets coefficients to 0)	arbitrarily picks 1 of a correlated group	lambda
Elastic net	\|\|y - X beta\|\|^2 + lambda * [(1 - alpha) * \|\|beta\|\|^2 / 2 + alpha * \|\|beta\|\|_1]	no (coordinate descent)	yes	keeps correlated groups together	lambda, alpha

Bias-variance trade-off in regularization

OLS is unbiased; any non-zero $\lambda$ introduces bias by pulling coefficients toward zero. The win is variance reduction: the shrinkage factor $s_j^{2}/(s_j^{2}+\lambda)$ is strictly below 1 on every direction, and tiny on small-singular-value directions. The optimal $\lambda$ is the one that balances the squared bias against the variance, and it is almost always strictly positive in finite samples [1].

Why a biased estimator can win

Gauss-Markov says OLS is best among unbiased linear estimators. Ridge is biased, so the theorem does not apply, and ridge often has lower MSE precisely because the bias-variance trade-off pays off when $X^{\top}X$ has small eigenvalues. This is the cleanest example in classical statistics of “controlled bias is good”.

Algorithm

Step	OLS via normal equations	OLS via gradient descent
1	Stack X (with intercept column) and y.	Stack X (with intercept column) and y.
2	Compute G = X^T X.	Initialise beta = 0.
3	Compute c = X^T y.	Choose step size eta below 2 / lambda_max(X^T X).
4	Solve G beta = c (e.g. Cholesky if G is PD).	Repeat: beta <- beta - eta * X^T (X beta - y) until \|\|grad\|\| is small.
5	Return beta_hat.	Return beta.
Cost	O(n d^2 + d^3) one-shot	O(n d) per step, no matrix inverse

# OLS via normal equations (pseudocode)
G = X.T @ X
c = X.T @ y
beta = solve(G, c)        # Cholesky / LU; do not call inv() in production
 
# OLS via gradient descent (pseudocode)
beta = zeros(d)
for step in range(max_steps):
    grad = X.T @ (X @ beta - y)
    beta = beta - eta * grad
    if norm(grad) < tol:
        break

Implementations

Three idiomatic implementations of the same model. The NumPy version solves the normal equations directly, the sklearn version is the production one-liner, and the PyTorch version uses autograd on the squared-error loss so the tensor view lines up with the math above.

python

import numpy as np

def ols_fit(X, y):
  """Closed-form OLS via the normal equations. Adds an intercept column."""
  n = X.shape[0]
  X1 = np.hstack([np.ones((n, 1)), X])
  G = X1.T @ X1
  c = X1.T @ y
  beta = np.linalg.solve(G, c)        # do not call inv() in production
  return beta[0], beta[1:]            # intercept, slopes

def ridge_fit(X, y, lam=1.0):
  """Closed-form ridge via (X^T X + lambda I)^{-1} X^T y. Penalty on slopes only."""
  n, d = X.shape
  X1 = np.hstack([np.ones((n, 1)), X])
  A = X1.T @ X1
  A[1:, 1:] += lam * np.eye(d)        # do not regularize the intercept
  beta = np.linalg.solve(A, X1.T @ y)
  return beta[0], beta[1:]

def ols_gd_fit(X, y, eta=1e-3, steps=2000):
  """OLS via batch gradient descent (no closed form)."""
  n, d = X.shape
  X1 = np.hstack([np.ones((n, 1)), X])
  beta = np.zeros(d + 1)
  for _ in range(steps):
      beta -= eta * (X1.T @ (X1 @ beta - y))
  return beta[0], beta[1:]

# small synthetic example
rng = np.random.default_rng(0)
X = rng.normal(size=(200, 3))
beta_true = np.array([1.5, -2.0, 0.5])
y = 0.7 + X @ beta_true + 0.3 * rng.normal(size=200)

print("OLS  :", ols_fit(X, y))
print("Ridge:", ridge_fit(X, y, lam=1.0))
print("GD   :", ols_gd_fit(X, y))

From-scratch reference. solve() is preferred over inv() for numerical stability.

Run the from-scratch NumPy normal-equations solver live in the browser:

Idle (packages: numpy)

(no output yet)

Real-world uses

Domain	System	Why linear regression
Real estate	Zillow Zestimate-style price models	Interpretable per-feature coefficients (sqft, beds, neighborhood) make legal and explainability requirements easier to meet.
Experimentation	A/B test impact estimation with covariates (CUPED-style)	Closed-form variance estimator and well-understood standard errors; covariate adjustment shrinks the variance of the treatment effect.
Finance	Fama-French and other factor models	Linear in factor exposures; OLS gives BLUE estimates of betas under the Gauss-Markov assumptions.
Pharma / clinical	Dose-response models	Regulatory pipelines require model forms whose assumptions can be checked, audited, and reproduced.
Causal inference	Difference-in-differences, instrumental variables, fixed effects	Most causal estimators are linear regressions in disguise, with closed-form standard errors.

Pitfalls

Multicollinearity inflates variance

If two columns of $X$ are nearly collinear, $X^{\top}X$ has a near-zero eigenvalue, and the OLS coefficient variance blows up along that direction. The point estimates are still unbiased but useless for inference. Diagnose with VIF; fix with ridge or by dropping a redundant feature.

Heteroskedasticity breaks standard errors

OLS point estimates remain unbiased under heteroskedasticity, but the textbook standard errors are wrong, so confidence intervals and p-values lie. Use heteroskedasticity-robust (sandwich, HC1, HC3) standard errors, or weighted least squares if you can model the variance structure.

Outliers dominate squared error

Squared error penalizes a residual of size $r$ by $r^{2}$ , so a single bad point can swing $\hat{\beta}$ noticeably. Use Huber regression for moderate contamination, RANSAC for adversarial outliers, or an explicit robust loss.

Standardize before regularizing

Ridge, lasso, and elastic net are not scale-invariant: a feature in millimetres gets penalized far more than the same feature in kilometres. Always standardize features (zero mean, unit variance) before running a regularized model, and never regularize the intercept.

Interview Q&A

Derive the OLS normal equations from the squared-error objective.research

Stack $n$ examples in $X \in \mathbb{R}^{n \times d}$ and $y \in \mathbb{R}^{n}$ , and define

\mathcal{L}(\beta) = \tfrac{1}{2}\,\lVert y - X\beta \rVert_{2}^{2} = \tfrac{1}{2}\,(y - X\beta)^{\top}(y - X\beta).

Expanding,

\mathcal{L}(\beta) = \tfrac{1}{2}\, y^{\top} y - \beta^{\top} X^{\top} y + \tfrac{1}{2}\,\beta^{\top} X^{\top} X \beta.

Differentiate with respect to $\beta$ using $\nabla_{\beta}(\beta^{\top} A \beta) = (A + A^{\top})\beta$ and $\nabla_{\beta}(c^{\top}\beta) = c$ :

\nabla_{\beta}\mathcal{L}(\beta) = -\,X^{\top} y + X^{\top} X\beta.

Set $\nabla_{\beta}\mathcal{L}(\beta) = 0$ to obtain $X^{\top} X \hat{\beta} = X^{\top} y$ , and if $X$ has full column rank, $\hat{\beta} = (X^{\top}X)^{-1} X^{\top} y$ . The Hessian is $X^{\top}X \succeq 0$ , so $\mathcal{L}$ is convex and the stationary point is the global minimum. Geometrically, $\hat{y} = X\hat{\beta} = P_X y$ is the orthogonal projection of $y$ onto $\operatorname{col}(X)$ , and the residual $r = y - \hat{y}$ satisfies $X^{\top} r = 0$ - which is the normal equations [1].

State and prove the Gauss-Markov theorem.research

Statement. If $y = X\beta + \varepsilon$ with $X$ non-random and full column rank, $\mathbb{E}[\varepsilon] = 0$ , and $\operatorname{Cov}(\varepsilon) = \sigma^{2} I$ , then among all linear unbiased estimators of $\beta$ , the OLS estimator $\hat{\beta} = (X^{\top}X)^{-1}X^{\top}y$ has the smallest variance (it is BLUE).

Proof sketch. Any linear estimator can be written $\tilde{\beta} = C y$ for a fixed $C \in \mathbb{R}^{d \times n}$ . Decompose $C = (X^{\top}X)^{-1}X^{\top} + D$ .

Unbiasedness. $\mathbb{E}[\tilde{\beta}] = C X \beta$ , so unbiasedness for all $\beta$ requires $C X = I$ . Substituting, $D X = 0$ .

Covariance. Using $\operatorname{Cov}(\varepsilon) = \sigma^{2} I$ ,

\operatorname{Cov}(\tilde{\beta}) = \sigma^{2}\, C C^{\top} = \sigma^{2}\bigl[(X^{\top}X)^{-1} + (X^{\top}X)^{-1}X^{\top}D^{\top} + D X (X^{\top}X)^{-1} + D D^{\top}\bigr].

The cross terms vanish because $D X = 0$ , leaving $\sigma^{2}\bigl[(X^{\top}X)^{-1} + D D^{\top}\bigr]$ . Since $D D^{\top} \succeq 0$ ,

\operatorname{Cov}(\tilde{\beta}) - \operatorname{Cov}(\hat{\beta}_{\mathrm{OLS}}) = \sigma^{2}\, D D^{\top} \succeq 0,

so $\hat{\beta}_{\mathrm{OLS}}$ has the smallest variance in the Loewner order. Equality requires $D = 0$ . ∎

The theorem says nothing about biased estimators, which is why ridge can have strictly lower MSE than OLS in finite samples [1].

Why is minimizing MSE equivalent to MLE under Gaussian noise?onsite

Assume $y_i = x_i^{\top}\beta + \varepsilon_i$ with $\varepsilon_i \sim \mathcal{N}(0, \sigma^2)$ i.i.d. Then

\log p(y \mid X, \beta, \sigma^2) = \sum_{i=1}^{n} \log \mathcal{N}(y_i \mid x_i^{\top}\beta, \sigma^2) = -\frac{n}{2}\log(2\pi\sigma^2) - \frac{1}{2\sigma^2}\sum_{i=1}^{n}(y_i - x_i^{\top}\beta)^2.

The first term does not depend on $\beta$ , and the second is a strictly decreasing function of the squared residual norm. Therefore

\arg\max_{\beta} \log p(y \mid X, \beta, \sigma^2) = \arg\min_{\beta}\,\lVert y - X\beta \rVert_{2}^{2}.

OLS is the MLE. The same calculation explains why ridge corresponds to MAP estimation under a Gaussian prior on $\beta$ , and why Laplace priors give lasso [2].

Derive the ridge closed form and explain why it works when X is rank-deficient.research

The ridge objective is

\mathcal{L}(\beta) = \tfrac{1}{2}\,\lVert y - X\beta \rVert_{2}^{2} + \tfrac{\lambda}{2}\,\lVert \beta \rVert_{2}^{2}.

Take the gradient and set to zero:

\nabla_{\beta}\mathcal{L}(\beta) = -X^{\top}y + X^{\top}X\beta + \lambda \beta = 0 \implies (X^{\top}X + \lambda I)\hat{\beta}_{\mathrm{ridge}} = X^{\top} y.

So $\hat{\beta}_{\mathrm{ridge}} = (X^{\top}X + \lambda I)^{-1} X^{\top} y$ .

Even if $X^{\top}X$ is singular (rank-deficient $X$ , or $d > n$ ), the matrix $X^{\top}X + \lambda I$ has all eigenvalues at least $\lambda > 0$ , so it is invertible and the solution is unique. This is why ridge is the standard fix for multicollinearity and high-dimensional underdetermined regression. The SVD view $X = U\Sigma V^{\top}$ makes the shrinkage explicit: ridge replaces $1/s_j$ with $s_j / (s_j^{2} + \lambda)$ , attenuating the small-singular-value directions where OLS variance is largest [4][1].

Check yourself

Score: 0 / 5(0 answered)

1. Which expression is the OLS estimator?
2. Under which assumption(s) is OLS the BLUE (Gauss-Markov)?
3. Why does ridge regression always have a unique closed-form solution?
4. Which property does lasso have that ridge does not?
5. What does it mean that OLS minimizing MSE equals MLE under Gaussian noise?

Build it yourself

A short evening project:

Generate synthetic data with a known $\beta^{\star}$ and Gaussian noise.
Fit OLS, ridge, lasso, and elastic net on the same training set across a grid of $\lambda$ values.
Plot validation MSE vs $\lambda$ and report the optimum for each model.
Add a few highly correlated features and re-run; watch ridge keep them all and lasso pick one.
Stretch: add deliberate heteroskedasticity, then refit OLS with HC3 robust standard errors and confirm the confidence intervals widen on high-variance regions.

Interactive Demo

Regression Fit Demo

Fitted line

y = 1.652x + 1.447

R² = 0.9005

What to try: Drag a point far from the cluster to see how one outlier pulls the regression line.

Regression line

Residuals

References

Hastie, T., Tibshirani, R., and Friedman, J. (2009). The Elements of Statistical Learning, 2nd ed.. Springer. https://hastie.su.domains/ElemStatLearn/
Bishop, C. M. (2006). Pattern Recognition and Machine Learning. Springer. ISBN 978-0-387-31073-2
Tibshirani, R. (1996). Regression Shrinkage and Selection via the Lasso. Journal of the Royal Statistical Society, Series B 58(1), 267-288. doi:10.1111/j.2517-6161.1996.tb02080.x
Hoerl, A. E. and Kennard, R. W. (1970). Ridge Regression: Biased Estimation for Nonorthogonal Problems. Technometrics 12(1), 55-67. doi:10.1080/00401706.1970.10488634
Zou, H. and Hastie, T. (2005). Regularization and variable selection via the elastic net. Journal of the Royal Statistical Society, Series B 67(2), 301-320. doi:10.1111/j.1467-9868.2005.00503.x

Perceptron Logistic Regression