From Least Squares to Likelihood: The Dual Nature of Linear Regression

Introduction

Imagine you want to predict the price $y$ of a house in City A. Naturally, the price depends on various factors — the size of the house, its location (proximity to schools, hospitals, etc.), its age, and perhaps a dozen more. In machine learning, we represent these influencing factors as a feature vector:$X = [x_1, x_2, ..., x_n]$, where each $x_i$ is a numerical representation of some property of the house.

Our goal is to learn a relationship between $X$ and $y$, so we can predict the price of new houses just from their features. One of the most classic ways to approach this is through linear regression:

This is the least squares perspective: find the weights $w_1, ..., w_n$ and bias $w_0$ that minimize the total squared error between predicted and actual prices across the training data.

In practice, no matter how many features we include, there will always be unpredictable or unmeasured influences — random fluctuations, measurement noise, and real-world complexities we choose to ignore for simplicity. That's why we don't write our model as $y = w^\top x$, but rather:

The $\varepsilon$ term captures this uncertainty — the deviation between the actual observed value and the model’s idealized prediction. It reminds us that the model is a simplification. As the statistician George E. P. Box famously said:

Thinking of $\varepsilon$ as a random variable — specifically, as a sample from a probability distribution — turns this from a deterministic equation into a probabilistic model. And as we’ll see next, this small conceptual shift leads us directly to maximum likelihood estimation.

The Least Squares Perspective

Before diving into the math, let’s address one notational point. In the previous section, we used $w_0$ to represent the intercept term. For clarity and to align with standard convention, we’ll now write the intercept as $b$, and keep the rest of the weights as a vector $w = [w_1, w_2, ..., w_n]^\top$. This allows us to express our model compactly as:

Now suppose we are given a dataset of $n$ examples, where each data point consists of a feature vector $x^{(i)} \in \mathbb{R}^d$ and a corresponding target value $y^{(i)} \in \mathbb{R}$. Our goal is to find the parameters $w$ and $b$ that best fit the data — in other words, that produce predictions $w^\top x^{(i)} + b$ that are as close as possible to the actual observed values $y^{(i)}$.

To make this concrete, here’s a snapshot of a real dataset used for regression: the California Housing dataset. Each row represents one example (a house), and each column is a numerical feature such as the median income of the area, average number of rooms, etc. The target value $y$ is the median house value.

In our notation, each row corresponds to a feature vector $x^{(i)}$, and the final column (the target) corresponds to $y^{(i)}$. If the dataset has $n$ rows and $d$ features, then:

• $$X \in \mathbb{R}^{n \times d}$$
  is the matrix of feature vectors
• $y \in \mathbb{R}^n$ is the vector of target values

To measure how good our predictions are, we define the squared error loss for each data point:

In machine learning, such loss functions are typically called objective functions — mathematical expressions we aim to minimize (or maximize) using optimization techniques like gradient descent or closed-form solvers. In this case, we use the squared error because it's:

• Differentiable, which makes it easy to optimize using calculus-based methods,
• Convex, so it guarantees a unique global minimum

Here's a visual intuition for what the squared error captures: Each blue point represents a house in our dataset, plotted against a single feature (average number of rooms). The red line is the regression model’s prediction. The dashed vertical lines show the errors — the distance between the true value and the predicted value for each data point.

The squared error loss function penalizes these vertical distances — especially large ones — and finds the line that minimizes the total penalty over the dataset.

The total error over the dataset is then the sum of these individual squared errors. This gives us the least squares objective function:

The optimal weights $w^*$ and intercept $b^*$ are those that minimize this loss:

This optimization problem defines the classical formulation of linear regression. It’s purely geometric and algebraic: we are finding the linear function that best fits the data by minimizing squared deviation.

To derive a closed-form solution, we express the entire dataset in matrix form. Let $X \in \mathbb{R}^{n \times d}$ be the matrix of feature vectors, where each row is a data point $x^{(i)\top}$ and $y \in \mathbb{R}^{n}$ be the vector of target values.

To incorporate the intercept $b$, we augment $X$with a column of ones, forming $\tilde{X} \in \mathbb{R}^{n \times (d+1)}$:

Similarly, we define the augmented parameter vector:

Now the prediction for all data points becomes simply:

And the least squares loss can be written compactly as:

This is a standard quadratic minimization problem. Let’s expand the squared loss function$\mathcal{L}(\tilde{w}) = \|y - \tilde{X} \tilde{w}\|_2^2$ into a fully expanded quadratic form:

Expanding this using standard matrix identities:

Before we compute the gradient, let's confirm that $\tilde{w}^\top \tilde{X}^\top \tilde{X} \tilde{w}$ is indeed a scalar (a single number), which justifies the symmetry trick used later in the derivation.

Recall the matrix dimensions:

• $\tilde{X} \in \mathbb{R}^{n \times (d+1)}$ — data matrix with bias column added
• $\tilde{w} \in \mathbb{R}^{(d+1) \times 1}$ — parameter vector

Then we have:

Here's the reasoning:

1. $\tilde{X} \tilde{w} \in \mathbb{R}^{n \times 1}$ — the predicted values
2. $\tilde{X}^\top (\tilde{X} \tilde{w}) \in \mathbb{R}^{(d+1) \times 1}$
3. $\tilde{w}^\top (\tilde{X}^\top \tilde{X} \tilde{w}) \in \mathbb{R}^{1 \times 1}$

Therefore, the whole expression collapses to a single real number. Since it's a scalar, we can freely take its transpose:

Which justifies the simplification in the gradient derivation.

We now take the gradient of $\mathcal{L}(\tilde{w})$ with respect to $\tilde{w}$. Using standard vector calculus identities (listed below), we obtain:

Setting the gradient to zero gives the first-order condition:

Cancelling the constant and rearranging, we get the normal equation:

And assuming $\tilde{X}^\top \tilde{X}$ is invertible, the unique solution is:

This is the closed-form solution to the least squares problem, derived directly from calculus and matrix algebra.

This is the closed-form solution to the linear regression problem. It gives us the optimal weights and intercept that minimize the total squared error across the training data.

But there’s also a powerful geometric interpretation here. If we look closely, we can see that linear regression is performing an orthogonal projection of the target vector $y \in \mathbb{R}^n$onto the column space of the design matrix $\tilde{X}$.

In other words, the predicted values $\hat{y} = \tilde{X} \tilde{w}$ live in the subspace spanned by the feature vectors (the columns of $\tilde{X}$). Among all possible vectors in this subspace, $\hat{y}$ is the one closest to the true output $y$ in Euclidean distance.

This makes $\hat{y}$ the best linear approximation of $y$ using the available features. The residual vector $y - \hat{y}$ is orthogonal to the subspace spanned by the feature vectors.