Normal Equations Derivation

With the matrix formulation established, we now face the central optimization problem: find the coefficient vector $\hat{\boldsymbol{\beta}}$ that minimizes the sum of squared errors.

Unlike many optimization problems in machine learning that require iterative algorithms, linear regression admits an exact, closed-form solution. This solution, derived from the normal equations, is one of the most elegant results in statistical learning theory.

The derivation we present here uses matrix calculus—specifically, derivatives of scalar functions with respect to vectors. This approach is not merely an exercise in mathematical sophistication; it's the foundation for understanding optimization throughout machine learning.

Recall from the previous page that we seek to minimize the sum of squared errors:

$$\hat{\boldsymbol{\beta}} = \underset{\boldsymbol{\beta}}{\arg\min} ; \text{SSE}(\boldsymbol{\beta}) = \underset{\boldsymbol{\beta}}{\arg\min} ; |\mathbf{y} - \mathbf{X}\boldsymbol{\beta}|^2$$

Expanding the squared norm:

$$\text{SSE}(\boldsymbol{\beta}) = (\mathbf{y} - \mathbf{X}\boldsymbol{\beta})^\top(\mathbf{y} - \mathbf{X}\boldsymbol{\beta})$$

This is a quadratic function of $\boldsymbol{\beta}$. Quadratic functions with positive semi-definite Hessian matrices have a unique global minimum (or a set of minima if the Hessian is singular). Our task is to find this minimum analytically.

Before deriving the normal equations, we need key results from matrix calculus. Let $\mathbf{x}$ be a column vector and $f(\mathbf{x})$ a scalar function. The gradient $ abla_{\mathbf{x}} f$ is a vector of partial derivatives:

$$ abla_{\mathbf{x}} f = \begin{bmatrix} \frac{\partial f}{\partial x_1} \ \frac{\partial f}{\partial x_2} \ \vdots \ \frac{\partial f}{\partial x_n} \end{bmatrix}$$

Here are the essential identities we'll use:

Proof sketch for the quadratic form gradient:

For $f(\mathbf{x}) = \mathbf{x}^\top\mathbf{A}\mathbf{x}$, we can expand: $$f(\mathbf{x}) = \sum_{i}\sum_{j} a_{ij} x_i x_j$$

Taking the partial derivative with respect to $x_k$: $$\frac{\partial f}{\partial x_k} = \sum_{j} a_{kj} x_j + \sum_{i} a_{ik} x_i = (\mathbf{A}\mathbf{x})_k + (\mathbf{A}^\top\mathbf{x})_k$$

Collecting terms: $ abla_\mathbf{x} f = (\mathbf{A} + \mathbf{A}^\top)\mathbf{x}$. When $\mathbf{A} = \mathbf{A}^\top$, this simplifies to $2\mathbf{A}\mathbf{x}$.

Let's expand the SSE objective to identify its structure:

$$\begin{aligned} \text{SSE}(\boldsymbol{\beta}) &= (\mathbf{y} - \mathbf{X}\boldsymbol{\beta})^\top(\mathbf{y} - \mathbf{X}\boldsymbol{\beta}) &= \mathbf{y}^\top\mathbf{y} - \mathbf{y}^\top\mathbf{X}\boldsymbol{\beta} - \boldsymbol{\beta}^\top\mathbf{X}^\top\mathbf{y} + \boldsymbol{\beta}^\top\mathbf{X}^\top\mathbf{X}\boldsymbol{\beta} \end{aligned}$$

Now, $\mathbf{y}^\top\mathbf{X}\boldsymbol{\beta}$ is a scalar. The transpose of a scalar equals itself, so: $$\mathbf{y}^\top\mathbf{X}\boldsymbol{\beta} = (\mathbf{y}^\top\mathbf{X}\boldsymbol{\beta})^\top = \boldsymbol{\beta}^\top\mathbf{X}^\top\mathbf{y}$$

Therefore, we can combine the middle terms: $$\text{SSE}(\boldsymbol{\beta}) = \mathbf{y}^\top\mathbf{y} - 2\boldsymbol{\beta}^\top\mathbf{X}^\top\mathbf{y} + \boldsymbol{\beta}^\top\mathbf{X}^\top\mathbf{X}\boldsymbol{\beta}$$

Identifying the components:

1. $\mathbf{y}^\top\mathbf{y}$ — A constant (doesn't depend on $\boldsymbol{\beta}$)
2. $-2\boldsymbol{\beta}^\top\mathbf{X}^\top\mathbf{y}$ — Linear in $\boldsymbol{\beta}$
3. $\boldsymbol{\beta}^\top\mathbf{X}^\top\mathbf{X}\boldsymbol{\beta}$ — Quadratic in $\boldsymbol{\beta}$

This is a quadratic function of $\boldsymbol{\beta}$ with:

• Linear coefficient: $-2\mathbf{X}^\top\mathbf{y}$
• Quadratic coefficient: $\mathbf{X}^\top\mathbf{X}$ (the Gram matrix)

The Gram matrix $\mathbf{X}^\top\mathbf{X}$ is symmetric and positive semi-definite, ensuring the objective is convex (bowl-shaped upward).

To find the minimum, we take the gradient with respect to $\boldsymbol{\beta}$ and set it to zero.

$$ abla_{\boldsymbol{\beta}} \text{SSE} = abla_{\boldsymbol{\beta}} \left[ \mathbf{y}^\top\mathbf{y} - 2\boldsymbol{\beta}^\top\mathbf{X}^\top\mathbf{y} + \boldsymbol{\beta}^\top\mathbf{X}^\top\mathbf{X}\boldsymbol{\beta} \right]$$

Applying our matrix calculus identities term by term:

1. $ abla_{\boldsymbol{\beta}}(\mathbf{y}^\top\mathbf{y}) = \mathbf{0}$ (constant term)

2. $ abla_{\boldsymbol{\beta}}(-2\boldsymbol{\beta}^\top\mathbf{X}^\top\mathbf{y}) = -2\mathbf{X}^\top\mathbf{y}$ (linear term)

3. $ abla_{\boldsymbol{\beta}}(\boldsymbol{\beta}^\top\mathbf{X}^\top\mathbf{X}\boldsymbol{\beta}) = 2\mathbf{X}^\top\mathbf{X}\boldsymbol{\beta}$ (quadratic term, using symmetry of $\mathbf{X}^\top\mathbf{X}$)

Combining: $$ abla_{\boldsymbol{\beta}} \text{SSE} = -2\mathbf{X}^\top\mathbf{y} + 2\mathbf{X}^\top\mathbf{X}\boldsymbol{\beta}$$

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

$$ abla_{\boldsymbol{\beta}} \text{SSE} = \mathbf{0}$$

$$-2\mathbf{X}^\top\mathbf{y} + 2\mathbf{X}^\top\mathbf{X}\boldsymbol{\beta} = \mathbf{0}$$

Dividing by 2 and rearranging:

$$\boxed{\mathbf{X}^\top\mathbf{X}\boldsymbol{\beta} = \mathbf{X}^\top\mathbf{y}}$$

These are the normal equations.

Solving for $\hat{\boldsymbol{\beta}}$:

If $\mathbf{X}^\top\mathbf{X}$ is invertible (i.e., has full rank), we can multiply both sides by $(\mathbf{X}^\top\mathbf{X})^{-1}$:

$$\boldsymbol{\beta} = (\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top\mathbf{y}$$

This gives us the OLS estimator:

$$\boxed{\hat{\boldsymbol{\beta}} = (\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top\mathbf{y}}$$

This is one of the most important formulas in statistics and machine learning. It expresses the optimal coefficients directly in terms of the data matrices.

Setting the gradient to zero finds stationary points—which could be minima, maxima, or saddle points. We must verify we have a minimum by examining the Hessian matrix (the matrix of second derivatives).

$$\mathbf{H} = abla^2_{\boldsymbol{\beta}} \text{SSE} = abla_{\boldsymbol{\beta}} \left( 2\mathbf{X}^\top\mathbf{X}\boldsymbol{\beta} - 2\mathbf{X}^\top\mathbf{y} \right) = 2\mathbf{X}^\top\mathbf{X}$$

The Hessian is constant (independent of $\boldsymbol{\beta}$) and equals $2\mathbf{X}^\top\mathbf{X}$.

Second-order sufficiency conditions:

For a stationary point to be a minimum, the Hessian must be positive semi-definite (for weak minimum) or positive definite (for strict, unique minimum).

Recall that $\mathbf{X}^\top\mathbf{X}$ is always positive semi-definite: $$\mathbf{v}^\top(\mathbf{X}^\top\mathbf{X})\mathbf{v} = (\mathbf{X}\mathbf{v})^\top(\mathbf{X}\mathbf{v}) = |\mathbf{X}\mathbf{v}|^2 \geq 0$$

Furthermore, $\mathbf{X}^\top\mathbf{X}$ is positive definite if and only if $\mathbf{X}$ has full column rank (i.e., no two columns are linearly dependent).

Conclusion: The stationary point $\hat{\boldsymbol{\beta}}$ is indeed a global minimum. If $\mathbf{X}^\top\mathbf{X}$ is invertible, it's the unique global minimum.

The OLS formula requires $(\mathbf{X}^\top\mathbf{X})^{-1}$ to exist, which happens when $\mathbf{X}^\top\mathbf{X}$ has full rank. This fails in several important scenarios:

Consequences:

When $\mathbf{X}^\top\mathbf{X}$ is singular:

1. Infinitely many solutions — The normal equations have a solution but it's not unique
2. Unstable estimates — Small data changes cause wild coefficient swings
3. Numerical issues — Even near-singularity causes precision problems

Solutions:

• Regularization: Ridge regression adds $\lambda\mathbf{I}$ to $\mathbf{X}^\top\mathbf{X}$, ensuring invertibility
• Feature selection: Remove redundant or perfectly correlated features
• Dimensionality reduction: Apply PCA before regression
• Moore-Penrose pseudoinverse: Provides a minimum-norm solution even when singular