Least Squares: Normal Equations & Regression

The Problem

A\mathbf{x} = \mathbf{b}

may have no solution —

\mathbf{b}

may not lie in the column space of

A

. This is typical when the system is overdetermined: more equations than unknowns, with the equations imposing contradictory constraints.

When no exact solution exists, the goal shifts from solving

A\mathbf{x} = \mathbf{b}

to minimizing the error:

\hat{\mathbf{x}} = \arg\min_{\mathbf{x}} \|A\mathbf{x} - \mathbf{b}\|^2

The quantity

\|A\mathbf{x} - \mathbf{b}\|^2 = \sum_i (A\mathbf{x} - \mathbf{b})_i^2

is the sum of squared residuals. The vector

\hat{\mathbf{x}}

that minimizes this sum is the least-squares solution.

Least Squares Notation

Notation

Least Squares Notation

A hat that marks a best guess rather than a solution, a named equation whose name means "perpendicular" and not "usual", and a superscript plus that stands in for an inverse the matrix does not have. All catalogued among the linear algebra symbols.

The estimate hat

\hat{\mathbf{x}}

,

\operatorname{proj}_W

and the projection matrix

P

come from projection notation;

\|\cdot\|

from norm notation;

A^{T}

and

A^{-1}

from matrix operation and inverse notation;

A = QR

from Gram-Schmidt notation;

\operatorname{Col}

from fundamental-space notation.

\hat{\mathbf{x}} = \arg\min_{\mathbf{x}} \|A\mathbf{x} - \mathbf{b}\|

x-hat is the x that minimises the residual length

The hat announces that this is a best available answer, not a solution: when

A\mathbf{x} = \mathbf{b}

has no solution,

\hat{\mathbf{x}}

is what replaces it. The arg min is doing the same job it does elsewhere — returning the input that minimises, not the minimum itself.

CasesSquaring is a notational convenience, not a change of problem:

\|A\mathbf{x} - \mathbf{b}\|^2

has the same minimiser as

\|A\mathbf{x} - \mathbf{b}\|

because squaring is increasing on non-negative numbers — and the square is what makes the derivative tractable, which is where the method's name comes from.

Do not confuseAn exact solution.

A\hat{\mathbf{x}} = \mathbf{b}

is generally false — that is the entire premise — so writing

\hat{\mathbf{x}}

without the hat quietly asserts something the setup denies.

A^{T}A\hat{\mathbf{x}} = A^{T}\mathbf{b}

the normal equations

A named system, and the name is geometric: "normal" here means perpendicular, recording that the residual is orthogonal to the column space — the condition The Normal Equations below derives. Multiplying through by

A^{T}

is what converts an unsolvable system into a solvable one.

CasesWhen

A^{T}A

is invertible the equations close into

\hat{\mathbf{x}} = (A^{T}A)^{-1}A^{T}\mathbf{b}

; that inverse exists exactly when

A

has independent columns, so the formula silently assumes full column rank and should not be written otherwise.

Do not confuse"Normal" as ordinary, or as the normal distribution. Neither is meant — the word points at the right angle in the geometry, and the collision with the statistical sense is pure accident of vocabulary, not a shared idea.

A^{+}

A-plus — the pseudoinverse of A

The plus superscript joins the corner census as a stand-in for an inverse that does not exist:

\hat{\mathbf{x}} = A^{+}\mathbf{b}

mirrors

\mathbf{x} = A^{-1}\mathbf{b}

in shape, and the changed symbol is the only warning that

A

may be non-square or singular. Its construction belongs to SVD notation.

Also written

A^{\dagger}

— the dagger — is equally standard and dominant in physics and numerical texts; both spellings mean the Moore–Penrose pseudoinverse, and the choice is house style rather than a difference in object.

Do not confuseThe superscript

+

of set theory and elsewhere. Here it neither adds nor denotes a positive part; the corner is a label, the same grammar as

A^{T}

and

A^{-1}

, and only its position identifies the job.

All linear algebra symbols →↑ Orthogonality

The Geometric Interpretation

The set of all vectors

A\mathbf{x}

as

\mathbf{x}

ranges over

\mathbb{R}^n

is the column space of

A

. Minimizing

\|A\mathbf{x} - \mathbf{b}\|

means finding the point in the column space closest to

\mathbf{b}

. That closest point is the orthogonal projection

\hat{\mathbf{b}} = \text{proj}_{\text{Col}(A)}\mathbf{b}

.

The least-squares solution

\hat{\mathbf{x}}

satisfies

A\hat{\mathbf{x}} = \hat{\mathbf{b}}

— it produces the projection, not the original

\mathbf{b}

. The residual

\mathbf{r} = \mathbf{b} - A\hat{\mathbf{x}} = \mathbf{b} - \hat{\mathbf{b}}

is the component of

\mathbf{b}

orthogonal to the column space. It lies in

\text{Col}(A)^\perp = \text{Null}(A^T)

.

The orthogonality condition

A^T\mathbf{r} = \mathbf{0}

— the residual is perpendicular to every column of

A

— is the geometric content of the least-squares solution. It is this condition that leads to the normal equations.

The residual, perpendicular to the column space

The target lies off the column space, so no exact solution exists. The best available answer is its projection, and what is left over is the residual — which comes out perpendicular to the column space rather than merely small. Minimising the error and making it perpendicular are the same requirement. See it for an inconsistent system on the least squares visualizer.

Everything algebraic that follows is a way of computing this projection without having to draw it.

The Normal Equations

The orthogonality condition

A^T(\mathbf{b} - A\hat{\mathbf{x}}) = \mathbf{0}

rearranges to

Normal Equations

A^TA\hat{\mathbf{x}} = A^T\mathbf{b}

Normal equations →

These are the normal equations. They form a square

n \times n

system regardless of the shape of

A

.

The matrix

A^TA

is always symmetric and positive semi-definite. When

A

has full column rank (the columns are linearly independent),

A^TA

is positive definite and invertible, giving a unique least-squares solution:

Least Squares Solution

\hat{\mathbf{x}} = (A^TA)^{-1}A^T\mathbf{b}

Least squares solution →

When

A

does not have full column rank,

A^TA

is singular and the normal equations have infinitely many solutions. All produce the same projection

\hat{\mathbf{b}} = A\hat{\mathbf{x}}

, but the

\hat{\mathbf{x}}

vectors differ. The minimum-norm solution — the one with smallest

\|\hat{\mathbf{x}}\|

— is selected by the pseudoinverse.

The normal equations assembled

Multiplying through by the transpose turns an unsolvable system into a solvable one, and the reason is the perpendicularity above: demanding that the residual be orthogonal to every column is exactly what these equations say. Form them for your own data on the least squares visualizer.

These equations are the classical route and the least numerically stable one, which is why QR is preferred in practice.

Worked Example: Fitting a Line

Fit a line

y = c_0 + c_1 x

to the data points

(1, 2)

,

(2, 3)

,

(3, 6)

,

(4, 7)

.

The model

y_i = c_0 + c_1 x_i

for each data point gives the system

A\mathbf{c} = \mathbf{y}

with

A = \begin{pmatrix} 1 & 1 \\ 1 & 2 \\ 1 & 3 \\ 1 & 4 \end{pmatrix}, \quad \mathbf{y} = \begin{pmatrix} 2 \\ 3 \\ 6 \\ 7 \end{pmatrix}

Four equations in two unknowns — overdetermined. The normal equations:

A^TA = \begin{pmatrix} 4 & 10 \\ 10 & 30 \end{pmatrix}, \quad A^T\mathbf{y} = \begin{pmatrix} 18 \\ 53 \end{pmatrix}

Solving:

\det(A^TA) = 120 - 100 = 20

.

(A^TA)^{-1} = \frac{1}{20}\begin{pmatrix} 30 & -10 \\ -10 & 4 \end{pmatrix}

.

\hat{\mathbf{c}} = \frac{1}{20}\begin{pmatrix} 30 & -10 \\ -10 & 4 \end{pmatrix}\begin{pmatrix} 18 \\ 53 \end{pmatrix} = \frac{1}{20}\begin{pmatrix} 540 - 530 \\ -180 + 212 \end{pmatrix} = \frac{1}{20}\begin{pmatrix} 10 \\ 32 \end{pmatrix} = \begin{pmatrix} 0.5 \\ 1.6 \end{pmatrix}

The best-fit line is

y = 0.5 + 1.6x

. The residuals are

2 - 2.1 = -0.1

,

3 - 3.7 = -0.7

,

6 - 5.3 = 0.7

,

7 - 6.9 = 0.1

. Their sum of squares

0.01 + 0.49 + 0.49 + 0.01 = 1.0

is the minimum achievable error for any line through these data.

A best-fit line through scattered points

The points do not lie on any line, so the system was never going to be consistent. What the fit returns is the line whose vertical errors have the smallest total square — and fitting a line is nothing more than least squares with two unknowns. Move the points and watch the fit respond on the least squares visualizer.

Fitting a parabola in the next section changes only the columns, never the method.

Worked Example: Fitting a Parabola

Fit a parabola

y = c_0 + c_1 x + c_2 x^2

to the same data

(1, 2)

,

(2, 3)

,

(3, 6)

,

(4, 7)

.

The design matrix gains a third column:

A = \begin{pmatrix} 1 & 1 & 1 \\ 1 & 2 & 4 \\ 1 & 3 & 9 \\ 1 & 4 & 16 \end{pmatrix}

The normal equations

A^TA\hat{\mathbf{c}} = A^T\mathbf{y}

now form a

3 \times 3

system. The machinery is identical — only the model matrix changes. A higher-degree model provides a closer fit (the residual sum of squares cannot increase when the model gains flexibility), but it also risks fitting noise rather than signal.

The framework generalizes to any linear model:

y = c_0 f_0(x) + c_1 f_1(x) + \cdots + c_k f_k(x)

where the functions

f_i

are chosen in advance. Each choice produces a different design matrix

A

, and the normal equations produce the best coefficients in the least-squares sense.

Model	Number of parameters	Row of the design matrix A	Typical use
Line: y = c₀ + c₁ x	2	[ 1, x_i ]	linear trends in one variable
Parabola: y = c₀ + c₁ x + c₂ x²	3	[ 1, x_i, x_i² ]	smooth one-variable curves with curvature
Polynomial of degree k	k + 1	[ 1, x_i, x_i², …, x_i^k ]	higher-order curve fitting; flexible but prone to overfit
Generic linear model: y = Σ_j c_j f_j(x)	k + 1	[ f₀(x_i), f₁(x_i), …, f_k(x_i) ]	arbitrary basis functions (sinusoids, exponentials, splines)
Multiple linear regression: y = β₀ + Σ_j β_j x_j	k + 1	[ 1, x_1,i, x_2,i, …, x_k,i ]	several predictors per observation

The Projection Matrix

The projection of

\mathbf{b}

onto the column space is

\hat{\mathbf{b}} = P\mathbf{b}

where

P = A(A^TA)^{-1}A^T

When the columns of

A

are orthonormal (

A = Q

with

Q^TQ = I

), this simplifies to

P = QQ^T

.

The projection matrix is symmetric (

P^T = P

) and idempotent (

P^2 = P

). The complementary matrix

I - P

projects onto the orthogonal complement

\text{Col}(A)^\perp

and extracts the residual:

\mathbf{r} = (I - P)\mathbf{b}

.

The minimum squared error is

\|\mathbf{r}\|^2 = \|(I - P)\mathbf{b}\|^2 = \|\mathbf{b}\|^2 - \|P\mathbf{b}\|^2

by the Pythagorean theorem, since

P\mathbf{b}

and

(I - P)\mathbf{b}

are orthogonal.

The Pseudoinverse

The matrix

A^+ = (A^TA)^{-1}A^T

(when

A

has full column rank) is called the left pseudoinverse of

A

. The least-squares solution is

\hat{\mathbf{x}} = A^+\mathbf{b}

.

The pseudoinverse satisfies

A^+A = I_n

(it is a left inverse of

A

), but

AA^+ \neq I_m

in general — the product

AA^+

equals the projection matrix

P

.

When

A

does not have full column rank, the Moore-Penrose pseudoinverse

A^+

is defined through the singular value decomposition: if

A = U\Sigma V^T

, then

A^+ = V\Sigma^+ U^T

, where

\Sigma^+

inverts the nonzero singular values and transposes the shape. The Moore-Penrose pseudoinverse gives the minimum-norm least-squares solution — the

\hat{\mathbf{x}}

of smallest length among all minimizers of

\|A\mathbf{x} - \mathbf{b}\|

.

Least Squares via QR

The QR decomposition

A = QR

provides a numerically superior method for solving the normal equations.

Substituting

A = QR

into

A^TA\hat{\mathbf{x}} = A^T\mathbf{b}

gives

R^TQ^TQR\hat{\mathbf{x}} = R^TQ^T\mathbf{b}

. Since

Q^TQ = I

, this simplifies to

R^TR\hat{\mathbf{x}} = R^TQ^T\mathbf{b}

. Multiplying both sides by

(R^T)^{-1}

:

Least Squares via QR

R\hat{\mathbf{x}} = Q^T\mathbf{b}

Least squares via qr →

The right-hand side

Q^T\mathbf{b}

is a vector of

n

dot products. The left-hand side is an upper triangular system, solved by back substitution.

The critical advantage is numerical. Forming

A^TA

explicitly squares the condition number of

A

, amplifying rounding errors. The QR approach works with

Q

and

R

directly, preserving the original conditioning. This is why QR-based least squares is the standard algorithm in numerical software, from MATLAB to Python's NumPy.

Regression as Least Squares

The entire framework of linear regression is a least-squares problem in matrix form.

Simple linear regression fits

y = \beta_0 + \beta_1 x + \epsilon

to data

(x_i, y_i)

. The design matrix has a column of ones and a column of

x_i

values. The normal equations produce the slope

\hat{\beta}_1

and intercept

\hat{\beta}_0

that minimize

\sum(y_i - \beta_0 - \beta_1 x_i)^2

.

Multiple linear regression fits

\mathbf{y} = X\boldsymbol{\beta} + \boldsymbol{\epsilon}

where

X

is the design matrix with rows for observations and columns for predictors. The normal equations

X^TX\hat{\boldsymbol{\beta}} = X^T\mathbf{y}

give the coefficient estimates.

The projection matrix

P = X(X^TX)^{-1}X^T

is called the hat matrix in statistics because it puts the "hat" on

\mathbf{y}

:

\hat{\mathbf{y}} = P\mathbf{y}

. The residual vector

\mathbf{e} = \mathbf{y} - \hat{\mathbf{y}} = (I - P)\mathbf{y}

is the component of

\mathbf{y}

orthogonal to the column space of

X

, and

\|\mathbf{e}\|^2

is the residual sum of squares.

Linear-algebra object	Statistics / regression name	Symbol and role
Coefficient matrix	design matrix / predictor matrix	X — rows = observations, columns = features (with an intercept column of 1s)
Vector of unknowns	coefficient (parameter) vector	β — parameters to be estimated
Right-hand side	response vector	y — observed outcomes
Least-squares solution	estimated coefficients	β̂ = (XᵀX)⁻¹ Xᵀ y
Projection matrix P	hat matrix	P = X(XᵀX)⁻¹Xᵀ — puts the "hat" on y
Projection of b onto Col(A)	fitted values	ŷ = P y = X β̂
Residual vector	residual (errors)	e = y − ŷ = (I − P) y
Squared residual length	residual sum of squares (RSS)	‖e‖² = Σᵢ (y_i − ŷ_i)²
Fraction of length captured by projection (centered)	coefficient of determination	R² — fraction of variation in y explained by the model

The Minimum Error

The least-squares error

\|\mathbf{r}\| = \|\mathbf{b} - A\hat{\mathbf{x}}\|

is the distance from

\mathbf{b}

to the column space of

A

. It is the length of the orthogonal component of

\mathbf{b}

with respect to

\text{Col}(A)

.

By the Pythagorean theorem, since

\hat{\mathbf{b}} = A\hat{\mathbf{x}}

and

\mathbf{r} = \mathbf{b} - \hat{\mathbf{b}}

are perpendicular:

\|\mathbf{b}\|^2 = \|\hat{\mathbf{b}}\|^2 + \|\mathbf{r}\|^2

The error is what remains after the projection accounts for as much of

\mathbf{b}

as possible. The ratio

\|\hat{\mathbf{b}}\|^2 / \|\mathbf{b}\|^2

(computed with centered data in regression) is the coefficient of determination

R^2

— the fraction of the total variation explained by the model. An

R^2

close to

1

means the column space captures nearly all of

\mathbf{b}

; close to

0

means the model explains little.

The error is zero if and only if

\mathbf{b} \in \text{Col}(A)

— if and only if the original system

A\mathbf{x} = \mathbf{b}

has an exact solution. In that case, the least-squares solution is the exact solution, and the two problems coincide.

Summary: Methods for Computing the Least-Squares Solution

The least-squares solution x̂ can be computed by several distinct routes, each with its own cost and numerical-stability profile. The table below collects the main methods alongside their formulas, operation counts, and stability characteristics, providing a single reference for choosing among them.

The four methods below all return the same

\hat{\mathbf{x}}

in exact arithmetic, so the choice between them is not about correctness. It is about a single decision: whether the method forms

A^{\mathsf{T}}A

or arranges never to see it. The grouping is that decision, and it is also the ranking.

Download

PNGRaster · 2×PDFPrint quality

Orthogonality · least squares

Four ways to compute the same solution

All four minimise $\|A\mathbf{x} - \mathbf{b}\|^2$ and return the same $\hat{\mathbf{x}}$ in exact arithmetic. In floating point they do not, and the difference is entirely about whether the method forms $A^{\mathsf{T}}A$ or avoids it.

4methods

Forms

A^{\mathsf{T}}A

— correct, fragile1

1

Normal equations§ 3

\sim mn^2 + n^3/3

; needs full column rank

A^{\mathsf{T}}A\hat{\mathbf{x}} = A^{\mathsf{T}}\mathbf{b}

The derivation everyone learns and the computation nobody should run. Setting the gradient to zero gives these equations directly, and

A^{\mathsf{T}}A

is symmetric positive definite so Cholesky solves it — but the squared condition number means the answer can lose half its digits before Cholesky begins.

Avoids it — the practical routes3

2

QR factorization§ 8

\sim 2mn^2

via Householder

R\hat{\mathbf{x}} = Q^{\mathsf{T}}\mathbf{b}

Factor

A = QR

, multiply through by

Q^{\mathsf{T}}

, and back-substitute on the triangular system.

Q

is orthogonal so it preserves lengths and therefore the condition number — the workhorse in every numerical library, and the reason QR is taught alongside least squares.

3

Pseudoinverse via SVD§ 7

no rank condition at all

\hat{\mathbf{x}} = V\Sigma^{+}U^{\mathsf{T}}\mathbf{b}

The only method that works when

A

is rank-deficient, where the others have no unique solution to find. Among the infinitely many minimisers it returns the one of smallest norm — a choice made deliberately rather than by accident of the algorithm. Costs more than QR and is worth it when the rank is in doubt.

4

Iterative methods§ 3

one matrix–vector product per step

Krylov iteration on

A^{\mathsf{T}}A\hat{\mathbf{x}} = A^{\mathsf{T}}\mathbf{b}

Never forms

A^{\mathsf{T}}A

explicitly — only its action on a vector, which costs two multiplications by

A

. The right choice when

A

is large and sparse and a factorization would destroy the sparsity. Convergence depends on conditioning, so preconditioning is part of the method rather than an optimisation.

Forming

A^{\mathsf{T}}A

squares the condition number — a matrix with

\kappa = 10^4

becomes

10^8

, and roughly half the available digits are lost before the solve begins. Every method below except the first is a way of not doing that, which is why the textbook derivation and the library implementation take different routes to the same answer.

Four ways to compute the same solution·/linear-algebra/orthogonality/least-squaresLearn Math Class

Why forming

A^{\mathsf{T}}A

hurts is worth stating in numbers. The condition number of

A^{\mathsf{T}}A

is the square of the condition number of

A

, so a matrix with

\kappa(A) = 10^4

— unremarkable in a fitting problem — produces a system with

\kappa = 10^8

. In double precision that consumes roughly half the available digits before the solve has started, and no amount of care in the solver recovers them.

QR avoids this because

Q

is orthogonal: multiplying by it preserves lengths, so it cannot amplify error, and the triangular system inherits

A

's conditioning rather than its square. The SVD goes further and drops the rank requirement entirely, which is the case where the other three have no unique answer to compute — among infinitely many minimisers it returns the one of least norm. That is a deliberate choice, not a default, and it is why the pseudoinverse is the right tool when the rank is uncertain.

Least Squares FAQ

Does $\hat{\mathbf{x}}$ satisfy $A\hat{\mathbf{x}} = \mathbf{b}$ ?

+

Generally no, and that failure is the whole premise. Least squares exists precisely because the system has no exact solution, so

\hat{\mathbf{x}}

is the vector making the residual as small as possible rather than zero. Dropping the hat and writing

A\mathbf{x} = \mathbf{b}

quietly asserts the very thing the setup denies.Read more →

Why minimise the squared norm instead of the norm itself?

+

Convenience, not a change of problem. Squaring is increasing on non-negative numbers, so

\|A\mathbf{x} - \mathbf{b}\|^2

and

\|A\mathbf{x} - \mathbf{b}\|

are minimised by the same

\hat{\mathbf{x}}

. What the square buys is a tractable derivative, since the square root would otherwise sit in the way. It also supplies the method's name.Read more →

Why are they called the "normal" equations?

+

Because normal here means perpendicular, pointing at the right angle in the geometry: the residual is orthogonal to the column space at the minimiser. It has nothing to do with ordinary in the everyday sense, and nothing to do with the normal distribution. That second collision is pure accident of vocabulary.Read more →

Least Squares

The Best Approximate Solution

The Problem

Least Squares Notation

Least Squares Notation

The Geometric Interpretation

The Normal Equations

Worked Example: Fitting a Line

Worked Example: Fitting a Parabola

The Projection Matrix

The Pseudoinverse

Least Squares via QR

Regression as Least Squares

The Minimum Error

Summary: Methods for Computing the Least-Squares Solution

Four ways to compute the same solution

Least Squares FAQ

Does x^\hat{\mathbf{x}}x^ satisfy Ax^=bA\hat{\mathbf{x}} = \mathbf{b}Ax^=b?

Why minimise the squared norm instead of the norm itself?

Why are they called the "normal" equations?

Does $\hat{\mathbf{x}}$ satisfy $A\hat{\mathbf{x}} = \mathbf{b}$ ?