# Normal Equations

The linear system whose solution is the least-squares fit, obtained by requiring the residual to be orthogonal to the column space of the design matrix.

**Rests on:** Matrix Multiplication, Orthogonality and Projections, Vector Spaces and Rank

**Leads to:** Least Squares Regression, Linear Regression, Partial Least Squares

*Mathematical Foundations · Linear Algebra*


## Why learn this

Imagine you are trying to predict next month’s electricity bill based on last month’s usage. You have data: for ten months, you recorded how much power was used (in kilowatt-hours) and what the bill actually was. You suspect a simple relationship: the bill is roughly proportional to usage. You want to find the single best multiplier that connects usage to cost. The method known as the normal equations gives you a direct recipe to calculate that multiplier from your data table. They are worth learning because they provide an exact answer without needing to guess or iterate, which is helpful when your dataset is small enough to fit on one screen and you need a clear, transparent calculation rather than a black-box approximation.


## The intuition

Let’s work through a concrete example. Suppose you have only two data points: Month 1 had 100 kWh usage and a $50 bill. Month 2 had 200 kWh usage and a $100 bill. You want to fit the line `Cost = m * Usage`. 

If you pick `m = 0.4`, your predictions are $40 (for 100 kWh) and $80 (for 200 kWh). The errors (the difference between predicted and actual) are -$10 and -$20. You want the error to be as small as possible in total. Specifically, you minimize the sum of the squared errors: `(-10)^2 + (-20)^2 = 100 + 400 = 500`.

Now try `m = 0.6`. Predictions are $60 and $120. Errors are +$10 and +$20. Sum of squares: `100 + 400 = 500` again.

Try `m = 0.5`. Predictions are $50 and $100. Errors are $0 and $0. Sum of squares is 0. This is clearly better.

The normal equations tell you how to find that perfect `m` without guessing. The formula for these equations, when applied to finding a single multiplier like this, effectively balances the data so that the calculation accounts for every piece of information equally. They work by ensuring that the remaining error (the part of the bill your simple line couldn't explain, often called the residual) has no connection to the input data (usage). In geometric terms, if you draw a vector representing all your actual bills, and another vector representing all your predicted bills, the gap between them must stand at a right angle to the space where your predictions live. If the gap weren’t perpendicular, you could tilt your line slightly to reduce the error further. The normal equations calculate exactly that tilting direction so the gap becomes perpendicular.

A key part of this process is understanding projection. When we say the prediction is 'projected' onto the space of possible usage values, we mean we are finding the closest possible match within our model's capabilities to the real data. Think of it like shining a light straight down on the usage line; the shadow cast by the actual bill amounts onto that line is your prediction.


## An analogy

Think of holding a string attached to a nail on a wall. The nail represents your actual data point (the true bill). You pull the other end of the string toward a straight ruler laid out on the floor beneath the nail. The ruler represents all possible predictions your model can make (since it's a simple line, it's a 1D space). When you let go, the string tightens and becomes perpendicular to the ruler. The point where the string touches the ruler is your 'best fit' prediction. The length of the remaining slack in the string is your error. If the string were at an angle to the ruler, you could slide your hand along the ruler to make the slack shorter. But when it’s perpendicular, any move along the ruler would only increase the length of the slack. That perpendicular point is exactly what the normal equations find.


## Definition

The normal equations are a specific set of calculations, represented by the formula (X^T X) beta = X^T y, which solve for the best-fit parameters by forcing the prediction error to be orthogonal (at right angles) to the input data features.

A residual is the difference between an observed value and the value predicted by the model.

To project a value onto a space is to find its closest corresponding point within that space.


## How it connects

This builds directly on your understanding of Matrix Multiplication and Linear Transformations, as these tools allow you to calculate how inputs map to predictions. It also relates to Orthogonality and Projections, because the 'perpendicular' condition means you are projecting the actual outcomes onto the space of possible predictions.


## The mathematics

The normal equations come from linear algebra, the branch of mathematics that deals with tables of numbers (matrices) and lists of numbers (vectors). Here, we use matrix multiplication to combine multiple inputs at once. The problem we want to solve: we have data — a set of input numbers (like hours studied) and matching output numbers (like exam scores). We want to find a straight line (or a hyperplane, in higher dimensions) that predicts the output from the input as well as possible. 'Well' means minimizing the sum of squared errors — the sum of the squares of the vertical distances from each data point to the line. The solution is found by setting up a matrix equation: the design matrix (let's call it X) holds the input data (with a column of 1s for the intercept), and the target vector y holds the outputs. We want to find a vector of weights w such that X·w is as close as possible to y. The 'as close as possible' condition is that the error vector (y − X·w) is perpendicular to every column of X. Perpendicular here means the dot product is zero. That condition, written out, gives the normal equations: Xᵀ·X·w = Xᵀ·y. In words: multiply the transpose of X with X, then with w, and set it equal to the transpose of X times y. Solving for w gives the best-fitting weights. We will work through this with real numbers next.

For a design matrix X (n rows, p columns) and target vector y (n rows),
the least-squares solution w* satisfies:

Xᵀ·X·w = Xᵀ·y

If Xᵀ·X is invertible, the solution is:

w* = (Xᵀ·X)⁻¹·Xᵀ·y


## A worked example

```
We have three data points: hours studied (x) and exam score (y).
(1, 2), (2, 3), (3, 5).
We want to fit a line: y = w₁·x + w₂ (where w₂ is the intercept, bias).
Design matrix X (with a column of 1s for the intercept):
X = [[1, 1], [1, 2], [1, 3]]   (first column is 1, second column is x)
Target vector y = [2, 3, 5].

Step 1: Step 1: Compute Xᵀ (transpose of X).
Xᵀ = [[1, 1, 1], [1, 2, 3]]   # The transpose swaps rows and columns. X had 3 rows and 2 columns, so Xᵀ has 2 rows and 3 columns.
Step 2: Step 2: Compute Xᵀ·X (a 2×2 matrix).
Xᵀ·X = [[1·1+1·1+1·1, 1·1+1·2+1·3], [1·1+2·1+3·1, 1·1+2·2+3·3]] = [[3, 6], [6, 14]]   # Matrix multiplication: take each row of the first matrix and each column of the second, multiply corresponding entries, and sum. This gives a 2×2 matrix because Xᵀ is 2×3 and X is 3×2.
Step 3: Step 3: Compute Xᵀ·y (a 2×1 vector).
Xᵀ·y = [[1·2+1·3+1·5], [1·2+2·3+3·5]] = [[10], [23]]   # Again multiply and sum: first row gives 2+3+5=10, second row gives 2+6+15=23.
Step 4: Step 4: Set up the normal equations: Xᵀ·X·w = Xᵀ·y.
[[3, 6], [6, 14]] · [w₁, w₂]ᵀ = [10, 23]ᵀ   # This is a system of two linear equations:
3w₁ + 6w₂ = 10
6w₁ + 14w₂ = 23
Step 5: Step 5: Solve for w₁ and w₂ using elimination.
From the first equation: w₁ = (10 − 6w₂)/3.
Substitute into second: 6·((10 − 6w₂)/3) + 14w₂ = 23
Simplify: 20 − 12w₂ + 14w₂ = 23 → 2w₂ = 3 → w₂ = 1.5
Then w₁ = (10 − 6·1.5)/3 = (10 − 9)/3 = 1/3 ≈ 0.333   # We solved the system to get the best-fit weights.

Result: The best-fit line is y = 0.333·x + 1.5. This is the least-squares line. For x=1, the prediction is 1.833, actual is 2 (error 0.167, squared 0.028); for x=2, prediction 2.167, actual 3 (error 0.833, squared 0.694); for x=3, prediction 2.5, actual 5 (error 2.5, squared 6.25). Sum of squared errors = 6.972. No other line gives a smaller sum.
```


## From scratch, in Python

```python
# We will solve the normal equations without any special libraries.
# First, we need a function to multiply two matrices (or a matrix and a vector).
# A matrix is a list of lists, e.g., [[1, 2], [3, 4]] is a 2x2 matrix.
# A vector is a list, e.g., [1, 2, 3].

def mat_mul(A, B):
    # A and B are lists of lists, or one is a list (vector).
    # We'll assume B is a matrix or a vector.
    if isinstance(B[0], list):
        # B is a matrix
        result = []
        for i in range(len(A)):
            row = []
            for j in range(len(B[0])):
                total = 0
                for k in range(len(A[0])):
                    total += A[i][k] * B[k][j]
                row.append(total)
            result.append(row)
        return result
    else:
        # B is a vector
        result = []
        for i in range(len(A)):
            total = 0
            for j in range(len(A[0])):
                total += A[i][j] * B[j]
            result.append(total)
        return result

def transpose(A):
    rows = len(A)
    cols = len(A[0])
    return [[A[i][j] for i in range(rows)] for j in range(cols)]

def invert_2x2(M):
    # Only works for 2x2 matrices. For larger, we'd need Gaussian elimination.
    a, b = M[0][0], M[0][1]
    c, d = M[1][0], M[1][1]
    det = a*d - b*c
    return [[d/det, -b/det], [-c/det, a/det]]

# Our data
X = [[1, 1], [1, 2], [1, 3]]  # first column is 1, second is x
y = [2, 3, 5]

# Compute Xᵀ
Xt = transpose(X)

# Compute Xᵀ·X
XtX = mat_mul(Xt, X)
print("Xᵀ·X =", XtX)

# Compute Xᵀ·y
Xty = mat_mul(Xt, y)
print("Xᵀ·y =", Xty)

# Invert Xᵀ·X (2x2)
XtX_inv = invert_2x2(XtX)
print("Inverse of Xᵀ·X =", XtX_inv)

# Multiply to get w* = (Xᵀ·X)⁻¹ · (Xᵀ·y)
w = mat_mul(XtX_inv, Xty)
print("Weights =", w)

# Prediction function
def predict(x):
    return w[0]*1 + w[1]*x  # intercept is w[0], slope is w[1]

print("Predictions:")
for point in X:
    x_val = point[1]
    print(f"x={x_val}, predicted={predict(x_val)}, actual={y[X.index(point)]}")
```


## The same thing, with a framework

```python
# Using NumPy, the standard library for numerical computing.
import numpy as np

# Data
X = np.array([[1, 1], [1, 2], [1, 3]])  # design matrix with intercept column
y = np.array([2, 3, 5])

# Solve the normal equations: (Xᵀ·X)⁻¹ · Xᵀ·y
XtX = X.T @ X
Xty = X.T @ y
w = np.linalg.inv(XtX) @ Xty
print("Weights:", w)

# Alternatively, numpy's lstsq function does the same more stably.
w_lstsq, _, _, _ = np.linalg.lstsq(X, y, rcond=None)
print("lstsq Weights:", w_lstsq)
```


## What the code does

The from-scratch code defines two main functions: `mat_mul` for matrix multiplication (looping over rows and columns) and `transpose`. It then sets up the design matrix `X` (each row is a data point: a 1 for the intercept and the x value). It computes `Xᵀ·X` and `Xᵀ·y` by hand, inverts the 2x2 matrix using the formula (determinant), and multiplies to get the weights. This is exactly the arithmetic you did by hand. The framework version using NumPy does the same thing in a few lines: `X.T` is the transpose, `@` is matrix multiplication, and `np.linalg.inv` computes the inverse. The `lstsq` function is the preferred way in practice because it is more numerically stable (it does not square the condition number).


## Common misconceptions

**The normal equations always give the best fit** — The normal equations give the line that is closest to the data. That must be the best model.

The normal equations find the best fit *for a linear function and for the least-squares loss*. If the true relationship is curved (e.g., advertising spend vs sales), a straight line may be poor. Also, the least-squares loss is sensitive to outliers; a different loss (like absolute error) has no closed-form solution. In fact: The normal equations give the exact solution to a specific optimization problem: minimizing the sum of squared residuals for a linear model. That is often a good starting point, but the 'best' model depends on the problem and the loss function.

**Residual means error in a bad sense** — Residual is a mistake, something wrong with the model.

In statistics, residual is just the difference between the observed value and the predicted value. It is not inherently bad; it is a natural part of any imperfect model. In fact: Residual is the vertical distance from the data point to the line. It is a measure of how much the model misses, but it is a technical term, not a judgment.

**You must invert XᵀX** — To solve the normal equations, I have to compute the inverse of XᵀX.

Inverting a matrix is almost always the wrong way to solve a linear system in practice. It is numerically unstable and slower. In fact: Computing the inverse is one way, but it is usually replaced by Gaussian elimination, QR decomposition, or SVD, which are more stable and faster.

**The normal equations are for single-variable regression only** — They work only when predicting one thing from one other thing.

The matrix form handles any number of features. The design matrix X can have many columns, and the same equations solve for multiple weights. In fact: The normal equations are general: they solve for any number of inputs, as long as the model is linear in the parameters.

**The line must pass through at least one data point** — A best-fit line should touch some of the observed points.

The least-squares line minimizes the total squared distance, and it often does not pass through any actual data point. In a noisy dataset, the line is an average trend, not an interpolating curve. In fact: The line is a summary of the trend; it does not have to hit any particular point.


## In short

The normal equations are a way to find the best-fitting line (or hyperplane) for a set of data points by minimizing the sum of squared errors. They come from a geometric idea: the errors (residuals) must be perpendicular to the space of possible predictions. The formula is: Xᵀ·X·w = Xᵀ·y. This is the exact solution that iterative methods approximate when the data is too large. Remember: it solves a specific optimization problem, and in practice you use QR or SVD instead of inverting XᵀX.


## Where this is used

**Large Language Models (LLMs):** The normal equations appear in the final linear layer that maps internal representations to token probabilities. During training, the model's embeddings (vectors) are fitted to predict the next token; this is a classification problem, not simple regression, but the concept of projecting onto a subspace appears in the matrix multiplications that make up the transformer layers. Specifically, the attention mechanism computes a weighted sum of value vectors — that is a projection onto the space of inputs. The normal equations give the closed-form solution for linear regression, which is a building block for understanding how embeddings are learned. A developer working on an LLM rarely uses normal equations directly (the data is huge, so they use gradient descent), but they use the same linear algebra: matrix multiplication of a design matrix (the input tokens) with weight matrices.

**Small Models (on-device):** When a model must run on a phone or a smartwatch, it must be small and fast. A common technique is to train a larger model and then compress it, or to fit a simple linear model directly. For example, a TinyML keyword spotter might use logistic regression (a linear model) to detect a wake word. The normal equations cannot be used directly because the data is streaming and there are many features, but they are the exact solution that iterative methods (like stochastic gradient descent) approximate. A developer might use a closed-form solution for a small linear model on a microcontroller, where the design matrix fits in memory. They would write code that builds the matrix and solves the normal equations, perhaps optimized in C.

**Retrieval-Augmented Generation (RAG):** RAG systems retrieve relevant documents by comparing embeddings — vectors representing text. The similarity is often the dot product between a query vector and document vectors. The concept of projecting a vector onto a subspace appears here: the 'column space' of the document matrix is the space of possible query representations. The normal equations are not directly used, but the geometric idea of orthogonality underlies how embeddings are structured: semantically similar items have high dot product (not perpendicular), and the embedding space is learned to make relevant documents close. A developer building a RAG system uses vector databases that compute dot products, and they think in terms of orthogonal directions when choosing embedding models.

**Agentic Systems (planning and tool use):** An agent that plans and calls tools often uses a model to predict the next action. The policy may be a neural network, but the underlying optimization during training is the same: minimize a loss function by adjusting weights. The normal equations are the closed-form solution for linear regression, which is the foundation of many reinforcement learning algorithms. For example, a simple agent might use linear function approximation to estimate the value of a state; solving for the weights that minimize the squared Bellman error is exactly a least-squares problem, and the normal equations give the exact solution for small state spaces. A developer implementing such an agent on a small scale could directly use the normal equations.

**Everyday Applied Machine Learning:** In any regression task — predicting house prices, forecasting sales, fraud detection (logistic regression) — the normal equations are the direct solution for linear models when the dataset is small to medium (thousands of rows, hundreds of features). A data scientist might use `numpy.linalg.lstsq` or the `LinearRegression` class from scikit-learn, which internally solve the normal equations (or a more stable variant). For example, to predict house prices from size and location, they build a design matrix with a column of 1s and the features, and the library computes the optimal weights. The developer writes:

```python
from sklearn.linear_model import LinearRegression
model = LinearRegression()
model.fit(X, y)
print(model.coef_)
```

The `fit` method solves the normal equations (using QR or SVD for stability). What breaks when done badly: if the features are highly correlated, the matrix XᵀX becomes ill-conditioned (the condition number is squared), leading to unstable weights. That is why practitioners use QR or SVD instead of directly inverting XᵀX.


---

Read this on the web: https://basavapurushottam.com/learn-ai/normal-equations
