Recap from Linear Algebra
With maximum likelihood estimation and
hypothesis testing in hand, we are ready to apply these tools to the most fundamental
predictive model in statistics. Linear regression sits at the intersection of linear algebra and probability. The algebraic formulation leads
to least-squares, while the probabilistic formulation leads to MLE. The two turn out to be equivalent under Gaussian noise.
Given observed data points \(\{(\boldsymbol{x}_i, y_i)\}_{i=1}^{n}\) where \(\boldsymbol{x}_i \in \mathbb{R}^d\) and \(y_i \in \mathbb{R}\), we assume the
linear regression model:
\[
\boldsymbol{y} = X\boldsymbol{\beta} + \boldsymbol{\epsilon}
\]
where \(X \in \mathbb{R}^{n \times d}\) is the design matrix, \(\boldsymbol{\beta} \in \mathbb{R}^d\) is the
parameter (weight) vector, \(\boldsymbol{y} \in \mathbb{R}^n\) is the observation vector, and
\(\boldsymbol{\epsilon} = \boldsymbol{y} - X\boldsymbol{\beta}\) is the error (noise) vector. For an estimate
\(\hat{\boldsymbol{\beta}}\), the observable vector \(\boldsymbol{y} - X\hat{\boldsymbol{\beta}}\) is the residual vector.
Here \(d\) is the number of features and \(n\) is the number of data points. The least-squares solution \(\hat{\boldsymbol{\beta}}\) must satisfy the
normal equations:
\[
X^\top X\hat{\boldsymbol{\beta}} = X^\top \boldsymbol{y}.
\]
They were derived geometrically, via orthogonal projection of \(\boldsymbol{y}\) onto \(\operatorname{Col}(X)\).
The model is "linear" in the parameters \(\boldsymbol{\beta}\), not in the features \(\boldsymbol{x}\). We can apply any
nonlinear transformation to the features. In a scalar case, for example, \(y = \beta_0 + \beta_1 x^2 + \beta_2 \sin(2\pi x)\) is still a
linear model because \(y\) is a linear combination of the transformed features with respect to \(\boldsymbol{\beta}\). We now show that the
same normal equations arise naturally from a probabilistic argument.
Linear Regression: A Probabilistic Perspective
We now reinterpret the linear model above probabilistically. Assume that each observation \(y_i\) is the realized value of a random variable
\(Y_i = \boldsymbol{x}_i^\top \boldsymbol{\beta} + \epsilon_i\), where the \(\epsilon_i\) are i.i.d. Gaussian noise of variance
\(\sigma^2 \gt 0\):
\[
\epsilon_i \sim \mathcal{N}(0, \sigma^2).
\]
Since \(\mathbb{E}[\epsilon_i] = 0\), the conditional mean of \(Y_i\) is
\[
\mathbb{E}[Y_i \mid \boldsymbol{x}_i] = \mu_i = \boldsymbol{x}_i^\top \boldsymbol{\beta}
\]
where \(\boldsymbol{\beta} \in \mathbb{R}^d\) contains the unknown regression parameters. This means that, given \(\boldsymbol{x}_i\), each
response is an independent Gaussian random variable
\(Y_i \sim \mathcal{N}(\boldsymbol{x}_i^\top \boldsymbol{\beta}, \sigma^2)\).
The conditional p.d.f. of a single observation \(y_i\) is therefore
\[
p(y_i \mid \boldsymbol{x}_i, \boldsymbol{\beta}, \sigma^2)
= \frac{1}{\sigma\sqrt{2\pi}} \exp\!\left\{-\frac{1}{2\sigma^2}(y_i - \boldsymbol{x}_i^\top \boldsymbol{\beta})^2\right\}
\]
and, by independence, the likelihood function for \(\boldsymbol{\beta}\) (with \(\sigma^2\) fixed) is
\[
\begin{align*}
L(\boldsymbol{\beta})
&= \prod_{i=1}^n p(y_i \mid \boldsymbol{x}_i, \boldsymbol{\beta}, \sigma^2) \\\\
&= \left(\frac{1}{\sigma\sqrt{2\pi}}\right)^{\!n} \exp\!\left\{-\frac{1}{2\sigma^2}\sum_{i=1}^n (y_i - \boldsymbol{x}_i^\top \boldsymbol{\beta})^2\right\}.
\end{align*}
\]
The log-likelihood function is
\[
\begin{align*}
\ln L(\boldsymbol{\beta})
&= n \ln \left(\frac{1}{\sigma \sqrt{2\pi}} \right) - \frac{1}{2\sigma^2} \sum_{i = 1}^n (y_i - \boldsymbol{x}_i^\top \boldsymbol{\beta})^2 \\\\
&= -\frac{n}{2}\ln(2\pi\sigma^2) - \frac{1}{2\sigma^2}\sum_{i=1}^n (y_i - \boldsymbol{x}_i^\top \boldsymbol{\beta})^2.
\end{align*}
\]
We set the gradient of the log-likelihood with respect to \(\boldsymbol{\beta}\) equal to zero. By the chain rule, the derivative of
\((y_i - \boldsymbol{x}_i^\top \boldsymbol{\beta})^2\) with respect to \(\boldsymbol{\beta}\) is
\(-2(y_i - \boldsymbol{x}_i^\top \boldsymbol{\beta})\boldsymbol{x}_i\). Combined with the outer coefficient \(-\frac{1}{2\sigma^2}\) from the
log-likelihood, the gradient simplifies to the overall factor \(\frac{1}{\sigma^2}\) below:
\[
\nabla_{\boldsymbol{\beta}} \ln L(\boldsymbol{\beta}) = \frac{1}{\sigma^2} \sum_{i=1}^n (y_i - \boldsymbol{x}_i^\top \boldsymbol{\beta}) \boldsymbol{x}_i = \mathbf{0}.
\]
Since \(\sigma^2 \gt 0\), the equation reduces to
\(\sum_{i=1}^n (y_i - \boldsymbol{x}_i^\top \boldsymbol{\beta}) \boldsymbol{x}_i = \mathbf{0}\). Distributing the sum and using the
scalar-vector identity
\((\boldsymbol{x}_i^\top \boldsymbol{\beta}) \boldsymbol{x}_i = \boldsymbol{x}_i \boldsymbol{x}_i^\top \boldsymbol{\beta}\):
\[
\begin{align*}
\sum_{i=1}^n y_i \boldsymbol{x}_i
&= \sum_{i=1}^n (\boldsymbol{x}_i^\top \boldsymbol{\beta}) \boldsymbol{x}_i \\\\
&= \sum_{i=1}^n \boldsymbol{x}_i \boldsymbol{x}_i^\top \boldsymbol{\beta} \\\\
&= X^\top X \boldsymbol{\beta}.
\end{align*}
\]
Recognizing the left side as \(X^\top \boldsymbol{y}\), we obtain
\[
X^\top X \boldsymbol{\beta} = X^\top \boldsymbol{y}.
\]
We have recovered exactly the normal equations of the linear-algebra derivation. When \(X^\top X\) is invertible, the MLE solution for linear
regression under Gaussian noise coincides with the least-squares solution:
\[
\begin{align*}
\hat{\boldsymbol{\beta}}_{\text{MLE}} &= \Big(\sum_{i=1}^n \boldsymbol{x}_i\boldsymbol{x}_i^\top \Big)^{-1}\Big(\sum_{i=1}^n \boldsymbol{x}_i y_i \Big) \\\\
&= (X^\top X)^{-1}X^\top \boldsymbol{y} \\\\
&= \hat{\boldsymbol{\beta}}_{\text{LS}}.
\end{align*}
\]
Verification via Direct Minimization
Consider the linear model \(\boldsymbol{y} = X\boldsymbol{\beta} + \boldsymbol{\epsilon}\) where \(X \in \mathbb{R}^{n \times d}\), \(\boldsymbol{y} \in \mathbb{R}^n\), and
\(\boldsymbol{\beta} \in \mathbb{R}^d\).
To obtain the least-squares solution \(\hat{\boldsymbol{\beta}}_{\text{LS}}\), we minimize the least-squares error objective function \(f(\boldsymbol{\beta})\):
\[
\begin{align*}
\hat{\boldsymbol{\beta}}_{\text{LS}} &= \arg \min_{\boldsymbol{\beta}} \| \boldsymbol{y} - X \boldsymbol{\beta} \|_{2}^2 \\\\
&= \arg \min_{\boldsymbol{\beta}} (\boldsymbol{y} - X \boldsymbol{\beta})^\top (\boldsymbol{y} - X \boldsymbol{\beta}) \\\\
&= \arg \min_{\boldsymbol{\beta}} \left( \boldsymbol{y}^\top \boldsymbol{y} - \boldsymbol{y}^\top X \boldsymbol{\beta} - \boldsymbol{\beta}^\top X^\top \boldsymbol{y} + \boldsymbol{\beta}^\top X^\top X \boldsymbol{\beta} \right) \\\\
&= \arg \min_{\boldsymbol{\beta}} f(\boldsymbol{\beta}).
\end{align*}
\]
We take the gradient of \(f(\boldsymbol{\beta})\) with respect to \(\boldsymbol{\beta}\) term by term. Fix a vector
\(\boldsymbol{b} \in \mathbb{R}^d\) and a matrix \(A \in \mathbb{R}^{d \times d}\), and write
\(\boldsymbol{b}^\top \boldsymbol{\beta} = \sum_k b_k \beta_k\) and
\(\boldsymbol{\beta}^\top A \boldsymbol{\beta} = \sum_{k, l} \beta_k A_{kl} \beta_l\). Differentiating in \(\beta_k\) yields
\(\nabla_{\boldsymbol{\beta}}(\boldsymbol{b}^\top \boldsymbol{\beta}) = \boldsymbol{b}\). The same computation gives
\(\nabla_{\boldsymbol{\beta}}(\boldsymbol{\beta}^\top A \boldsymbol{\beta}) = (A + A^\top)\boldsymbol{\beta}\), because \(\beta_k\)
occurs in both factors of the quadratic form. The four terms of \(f\) then have the gradients:
\[
\begin{align*}
\nabla_{\boldsymbol{\beta}}(\boldsymbol{y}^\top \boldsymbol{y}) &= \mathbf{0}, \\\\
\nabla_{\boldsymbol{\beta}}(\boldsymbol{y}^\top X \boldsymbol{\beta}) &= X^\top \boldsymbol{y}, \\\\
\nabla_{\boldsymbol{\beta}}(\boldsymbol{\beta}^\top X^\top \boldsymbol{y}) &= X^\top \boldsymbol{y}, \\\\
\nabla_{\boldsymbol{\beta}}(\boldsymbol{\beta}^\top X^\top X \boldsymbol{\beta}) &= 2 X^\top X \boldsymbol{\beta} \quad \text{(since \(X^\top X\) is symmetric).}
\end{align*}
\]
Assembling the four terms and setting the gradient to zero:
\[
\begin{align*}
\nabla_{\boldsymbol{\beta}} f(\boldsymbol{\beta})
&= - X^\top \boldsymbol{y} - X^\top \boldsymbol{y} + 2X^\top X\boldsymbol{\beta} \\\\
&= -2 X^\top \boldsymbol{y} + 2X^\top X\boldsymbol{\beta} \\\\
&= \mathbf{0},
\end{align*}
\]
which simplifies to the normal equations \(X^\top X\boldsymbol{\beta} = X^\top \boldsymbol{y}\). Since
\(\boldsymbol{\beta}^\top X^\top X \boldsymbol{\beta} = \| X \boldsymbol{\beta} \|_{2}^2 \geq 0\), the matrix \(X^\top X\) is
positive semi-definite, so \(f\) is convex and every stationary point is a global minimizer. When \(X^\top X\) is invertible, the
unique minimizer is
\[
\hat{\boldsymbol{\beta}}_{\text{LS}} = (X^\top X)^{-1}X^\top \boldsymbol{y}.
\]
The matrix \(X^\top X\) is symmetric, and it is invertible if and only if \(\operatorname{rank}(X) = d\), that is, the columns of \(X\)
are linearly independent (this requires \(n \geq d\)).
Connections to Machine Learning
The equivalence \(\hat{\boldsymbol{\beta}}_{\text{MLE}} = \hat{\boldsymbol{\beta}}_{\text{LS}}\) reveals a fundamental principle.
Minimizing squared error is equivalent to maximum likelihood under Gaussian noise. The connection extends broadly.
Adding a squared \(\ell^2\) penalty
(ridge regression)
corresponds to placing a zero-mean Gaussian prior on \(\boldsymbol{\beta}\) and computing
the MAP estimate, while an \(\ell^1\) penalty
(Lasso) corresponds to a
zero-mean Laplace prior. The
bias-variance decomposition
of the ridge prediction error mirrors the split of the mean squared error
into variance and squared bias, which we met with maximum likelihood estimation.
When the Gaussian noise assumption fails, one may use generalized linear models or robust regression methods that replace the squared
loss with alternatives derived from other likelihood models.
Interactive Statistical Regression Tool
The tool below fits the model of this page, either to uploaded or pasted CSV data or to the built-in samples, and adds the standard
inference that the Gaussian noise model supports; those inferential formulas are stated here without derivation. The coefficients are
computed by the QR route. Under the Gaussian noise model above, with \(n \gt p\), the entries of the coefficient table are exact
finite-sample inference. The standard error of each coefficient is the square root of the corresponding diagonal entry of
\(\hat{\sigma}^2 (X^\top X)^{-1}\), with \(\hat{\sigma}^2 = \text{SSE}/(n - p)\) for \(p\) estimated coefficients, and the \(t\)-statistics
with their two-sided \(p\)-values and the 95% confidence intervals are built from the \(t_{n-p}\) distribution. The critical value is
computed for the actual sample size, not approximated by the large-sample 1.96.
The summary line reports \(R^2\), adjusted \(R^2\), and the overall \(F\)-test of the hypothesis that every non-intercept coefficient is
zero. That test applies the hypothesis-testing machinery to the model as a whole.
In single-predictor mode, the regression plot distinguishes two different 95% bands: the confidence band for the mean
response \(\mathbb{E}[Y \mid \boldsymbol{x}]\), and the wider prediction band for a new observation. The prediction band is
wider because its variance carries an extra \(\sigma^2\), the irreducible noise that no amount of data removes.
The diagnostics panel checks the Gaussian assumption that made MLE equal least squares in the first place. The residuals-vs-fitted plot
exposes nonlinearity and heteroscedasticity (try the "Heteroscedastic" sample and watch the funnel), and the normal QQ plot puts the
standardized residuals against Gaussian quantiles (try "Outliers" and watch the tails peel off the line). When predictors are exactly
collinear or \(n \lt p\), the matrix \(X^\top X\) is singular, and the tool says so instead of printing meaningless coefficients. The check
is a numerical version of the invertibility condition \(\operatorname{rank}(X) = d\) from the proof above: a column that is linearly
dependent on the earlier ones up to a small relative tolerance counts as dependent, so designs that are singular to within that tolerance
are declined as well. When \(n = p\) and \(X\) has full rank, \(X^\top X\) is invertible, but the fit interpolates the data exactly. The
estimate \(\hat{\sigma}^2 = \text{SSE}/(n - p)\) is then undefined, so the tool reports the coefficients without standard errors or
\(p\)-values. It does the same when \(n \gt p\) but the residuals vanish up to a small relative tolerance (as in a perfect fit, which
includes a constant response), since \(\hat{\sigma}^2\) is then negligible.