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 - and 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 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 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}.
\]
These were derived geometrically in
Least-Squares Problems
via orthogonal projection of \(\boldsymbol{y}\) onto \(\operatorname{Col}(X)\).
An important point: the model is "linear" in the parameters \(\boldsymbol{\beta}\), not in the features \(\boldsymbol{x}\).
We can apply any nonlinear transformation to the features - for example, in a scalar case,
\(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\) that depends on the predictor \(\boldsymbol{x}_i\), corrupted by i.i.d. Gaussian noise:
\[
\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] = \mu_i = \boldsymbol{x}_i^\top \boldsymbol{\beta}
\]
where \(\boldsymbol{\beta} \in \mathbb{R}^d\) contains the unknown regression parameters.
This means 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 > 0\), this 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}.
\]
This is exactly the normal equations from the linear-algebra derivation. Thus, 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*}
\]
Taking the gradient of \(f(\boldsymbol{\beta})\) with respect to \(\boldsymbol{\beta}\) term by term,
using the standard matrix-calculus identities \(\nabla_{\boldsymbol{\beta}}(\boldsymbol{b}^\top \boldsymbol{\beta}) = \boldsymbol{b}\)
and \(\nabla_{\boldsymbol{\beta}}(\boldsymbol{\beta}^\top A \boldsymbol{\beta}) = (A + A^\top)\boldsymbol{\beta}\):
\[
\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} \qquad \text{(since \(X^\top X\) is symmetric).}
\end{align*}
\]
Assembling the four terms and setting the gradient to zero:
\[
\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},
\]
which simplifies to the normal equations \(X^\top X\boldsymbol{\beta} = X^\top \boldsymbol{y}\).
When \(X^\top X\) is invertible, the unique minimizer is
\[
\hat{\boldsymbol{\beta}}_{\text{LS}} = (X^\top X)^{-1}X^\top \boldsymbol{y}.
\]
Note that \(X^\top X\) is symmetric and positive semi-definite. It is invertible if and only if \(\operatorname{rank}(X) = d\), i.e.,
the columns of \(X\) are linearly independent (which 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. This connection extends broadly.
Adding an \(\ell^2\) penalty (ridge regression)
corresponds to placing a Gaussian prior on \(\boldsymbol{\beta}\) and computing the MAP estimate, while
an \(\ell^1\) penalty (Lasso) corresponds to a
Laplace prior.
The bias-variance tradeoff
in regularized regression directly reflects the MSE decomposition developed
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 runs the full statistical analysis this page develops — on your own data (upload or
paste a CSV) or on the built-in samples. The coefficients are computed by the QR route from
the least-squares page, and every entry of the coefficient table is genuine inference: standard errors
from \(\hat{\sigma}^2 (X^\top X)^{-1}\) with \(\hat{\sigma}^2 = \text{SSE}/(n - p)\),
\(t\)-statistics with their two-sided \(p\)-values, and 95% confidence intervals built from the
\(t_{n-p}\) distribution — the critical value is computed for your 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, connecting directly to the
hypothesis-testing machinery of the previous page.
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,
whose extra width is exactly the irreducible noise \(\sigma^2\) 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 collinear or \(n \leq p\), the tool reports that \(X^\top X\) is singular instead of
printing meaningless coefficients — the invertibility condition \(\operatorname{rank}(X) = d\) from the
proof above, enforced live.