3.5.3 Least squares | SigmaField

A system with no solution

Lesson 3.3.2 solved a tall system, three equations in two unknowns. Change its last right-hand side from 4 to 5:

\begin{aligned} x_1 + x_2 &= 3 \\ x_1 - x_2 &= 1 \\ x_1 + 2x_2 &= 5 \end{aligned}

The first two equations give $(2, 1)$ , and then the third reads $4 = 5$ . As a matrix equation this is $T\mathbf{x} = \mathbf{b}$ with $\mathbf{b} = (3, 1, 5)$ . Lesson 3.5.1 found that $\mathbf{b}$ is not in the column space of $T$ , so no $\mathbf{x}$ solves it.

A system with more equations than unknowns is overdetermined. Data gives such systems all the time: many measurements, few unknowns, and a little error in each measurement.

The best answer

When no $\mathbf{x}$ makes $T\mathbf{x} = \mathbf{b}$ , look for the $\mathbf{x}$ that makes $T\mathbf{x}$ as close to $\mathbf{b}$ as possible. The gap $T\mathbf{x} - \mathbf{b}$ has one entry per equation, the amount by which that equation misses. Its squared length is the sum of the squared misses.

$\mathbf{x}$	$T\mathbf{x}$	$T\mathbf{x} - \mathbf{b}$	$\lVert T\mathbf{x} - \mathbf{b} \rVert^2$
$(2, 1)$	$(3, 1, 4)$	$(0, 0, -1)$	1
$(2, 1.5)$	$(3.5, 0.5, 5)$	$(0.5, -0.5, 0)$	0.5
$(2.5, 1)$	$(3.5, 1.5, 4.5)$	$(0.5, 0.5, -0.5)$	0.75
$(2.1429, 1.2857)$	$(3.4286, 0.8571, 4.7143)$	$(0.4286, -0.1429, -0.2857)$	0.2857

The first row meets two equations exactly and misses the third by 1. The last row meets none of them exactly, yet its total is the smallest. It is the least-squares solution, found below.

Rule: Least squares

\hat{\mathbf{x}} \text{ makes } \lVert A\mathbf{x} - \mathbf{b} \rVert^2 \text{ smallest}

$\hat{\mathbf{x}}$ , read "x hat", is the least-squares solution of $A\mathbf{x} = \mathbf{b}$ . The squared length adds the squares of the misses. When the system has a solution, $\hat{\mathbf{x}}$ is that solution, with every miss 0.

The normal equations

Every product $A\mathbf{x}$ lies in the column space of $A$ . The point of the column space closest to $\mathbf{b}$ is its projection $\mathbf{p}$ (lesson 3.5.2). So $A\hat{\mathbf{x}} = \mathbf{p}$ , and the residual $\mathbf{b} - A\hat{\mathbf{x}}$ is orthogonal to every column of $A$ .

Rule: The normal equations

\begin{gathered} A^\top(\mathbf{b} - A\hat{\mathbf{x}}) = \mathbf{0} \\[4pt] A^\top A\,\hat{\mathbf{x}} = A^\top\mathbf{b} \end{gathered}

"Normal" means at right angles. When the columns of $A$ are independent, $A^\top A$ has an inverse and $\hat{\mathbf{x}} = (A^\top A)^{-1}A^\top\mathbf{b}$ , exactly the projection formula of lesson 3.5.2.

For the system above:

\begin{gathered} T^\top T = \begin{pmatrix} 3 & 2 \\ 2 & 6 \end{pmatrix}, \quad T^\top\mathbf{b} = \begin{pmatrix} 9 \\ 12 \end{pmatrix} \\[6pt] 3x_1 + 2x_2 = 9, \quad 2x_1 + 6x_2 = 12 \end{gathered}

Three times the first equation minus the second gives $7x_1 = 15$ . So $\hat{\mathbf{x}} = \left(\frac{15}{7}, \frac{9}{7}\right) \approx (2.1429, 1.2857)$ .

The residual is $\mathbf{b} - T\hat{\mathbf{x}} = \frac{1}{7}(-3, 1, 2)$ . It is orthogonal to both columns: $\frac{1}{7}(-3 + 1 + 2) = 0$ and $\frac{1}{7}(-3 - 1 + 4) = 0$ . Its squared length is $\frac{14}{49} = \frac{2}{7} \approx 0.2857$ , the smallest in the table.

T

import numpy as np T = np.array([[1.0, 1.0], [1.0, -1.0], [1.0, 2.0]]) b = np.array([3.0, 1.0, 5.0]) print(T.T @ T, T.T @ b) x_hat = np.linalg.solve(T.T @ T, T.T @ b) r = b - T @ x_hat print(np.round(x_hat, 4), np.round(r, 4)) print(np.round(T.T @ r, 10) + 0.0, round(r @ r, 4)) for x in [(2, 1), (2, 1.5), (2.5, 1), tuple(x_hat)]: miss = T @ np.array(x) - b print(np.round(x, 4), np.round(miss, 4), round(miss @ miss, 4))

⌘/Ctrl + Enter runs the cell

Fitting a line to points

Lesson 1.2.1 drew the line through two points. Points that do not lie on one line have no line through them all, but least squares finds the line that comes closest.

Take five points: $(0, 1)$ , $(1, 2)$ , $(2, 4)$ , $(3, 4)$ and $(4, 6)$ . Look for a line $y = c + mx$ , with intercept $c$ and slope $m$ . Each point asks for $c + m x_i = y_i$ , so there are five equations in two unknowns:

\begin{pmatrix} 1 & 0 \\ 1 & 1 \\ 1 & 2 \\ 1 & 3 \\ 1 & 4 \end{pmatrix} \begin{pmatrix} c \\ m \end{pmatrix} = \begin{pmatrix} 1 \\ 2 \\ 4 \\ 4 \\ 6 \end{pmatrix}

The column of ones carries the intercept, and the column of $x$ values the slope. The normal equations hold the sums of the data:

\begin{gathered} A^\top A = \begin{pmatrix} 5 & 10 \\ 10 & 30 \end{pmatrix}, \quad A^\top\mathbf{y} = \begin{pmatrix} 17 \\ 46 \end{pmatrix} \end{gathered}

Here 5 is the number of points, 10 the sum of the $x_i$ and 30 the sum of their squares. On the right, 17 is the sum of the $y_i$ and 46 the sum of the products $x_iy_i$ .

Solving gives $c = 1$ and $m = 1.2$ : the line $y = 1 + 1.2x$ .

$x_i$	0	1	2	3	4
$y_i$	1	2	4	4	6
fitted $1 + 1.2x_i$	1	2.2	3.4	4.6	5.8
residual	0	$-0.2$	0.6	$-0.6$	0.2

The squared residuals add up to 0.8. No other line does better: the line through the first and last points leaves 0.875, and $y = 1 + x$ leaves 2.

Two lines through five points, with each residual drawn as a vertical gap. Least squares makes the squared gaps add up to the least possible, 0.8; the line through the end points leaves 0.875.

Each residual is a vertical gap, a miss in $y$ , not the distance from the point to the line at right angles. The right angle of lesson 3.5.2 lives in $\mathbb{R}^5$ : the residual vector is orthogonal to the two columns of $A$ .

That right angle says two things you can check. The residuals add up to 0, since they are orthogonal to the column of ones. And the line passes through the mean point $(\bar{x}, \bar{y}) = (2, 3.4)$ : $1 + 1.2 \times 2 = 3.4$ .

The slope as a projection

Centre the data as in case study 3.1.4: subtract the means to get $\tilde{\mathbf{x}}$ and $\tilde{\mathbf{y}}$ . The slope is the coefficient of the projection of $\tilde{\mathbf{y}}$ onto $\tilde{\mathbf{x}}$ (lesson 3.5.2), and the intercept puts the line through the mean point.

Rule: Least-squares line

\begin{gathered} m = \frac{\tilde{\mathbf{x}} \cdot \tilde{\mathbf{y}}}{\tilde{\mathbf{x}} \cdot \tilde{\mathbf{x}}} \\[4pt] c = \bar{y} - m\bar{x} \end{gathered}

$\tilde{\mathbf{x}} = \mathbf{x} - \bar{x}\,\mathbf{1}$ and $\tilde{\mathbf{y}} = \mathbf{y} - \bar{y}\,\mathbf{1}$ are the centred data. Here $\tilde{\mathbf{x}} \cdot \tilde{\mathbf{y}} = 12$ and $\tilde{\mathbf{x}} \cdot \tilde{\mathbf{x}} = 10$ , so $m = 1.2$ , and $c = 3.4 - 1.2 \times 2 = 1$ .

Case study 3.5.5 uses exactly this projection to find a hedge.

A

import numpy as np x = np.array([0, 1, 2, 3, 4], dtype=float) y = np.array([1, 2, 4, 4, 6], dtype=float) A = np.column_stack([np.ones(5), x]) print(A.T @ A, A.T @ y) c, m = np.linalg.solve(A.T @ A, A.T @ y) print(round(c, 4), round(m, 4)) res = y - (c + m * x) print(np.round(res, 4) + 0.0, round(res @ res, 4), round(res.sum(), 10) + 0.0, round(x @ res, 10) + 0.0) xc, yc = x - x.mean(), y - y.mean() print(xc @ yc, xc @ xc, xc @ yc / (xc @ xc)) print(np.round(np.polyfit(x, y, 1), 4))

⌘/Ctrl + Enter runs the cell

A noisy line

Real data has more points and more noise. The next cell makes 20 points from the line $y = 2 + 0.5x$ , adds a random error to each, and fits a line. The data is invented with a fixed seed, so every run gives the same points.

The fit is $y = 2.4162 + 0.4571x$ : close to the line the points came from, but not on it, because of the noise. How far off a fit like this can be is a question of statistics, which module 5.10 adds.

\hat{\mathbf{x}}

import numpy as np rng = np.random.default_rng(353) x = np.arange(1, 21, dtype=float) y = 2 + 0.5 * x + rng.normal(0, 1, 20) A = np.column_stack([np.ones(20), x]) x_hat, ss, rank, sv = np.linalg.lstsq(A, y, rcond=None) print(np.round(x_hat, 4), np.round(ss, 4), rank) print(np.round(np.linalg.solve(A.T @ A, A.T @ y), 4))

⌘/Ctrl + Enter runs the cell

Solve with lstsq, not with AᵀA

The normal equations are the idea; np.linalg.lstsq(A, b, rcond=None) is the tool. It works on $A$ itself and never forms $A^\top A$ .

Forming $A^\top A$ squares the condition number of lesson 3.3.2, which measures how much a system can magnify small errors such as rounding. A computer keeps about 16 significant digits, and a condition number of $10^k$ can cost about $k$ of them.

Calendar years. Fit a trend against the years 2015 to 2024. The column of ones and the column of years point almost the same way: the angle between them is 0.08°. Then $\operatorname{cond}(A) \approx 1.4 \times 10^6$ , but $\operatorname{cond}(A^\top A) \approx 2.0 \times 10^{12}$ , which leaves only about 4 safe digits.

Centring the years first, as for the slope above, makes the two columns orthogonal. The condition number drops to 2.8723. Centring helps lstsq too: practice set 3.5.4 meets $x$ values so large that it loses the slope without it.

Columns that differ in the eighth decimal. Take columns $(1, 10^{-8}, 0)$ and $(1, 0, 10^{-8})$ . Each diagonal entry of $A^\top A$ is $1 + 10^{-16}$ , which a computer rounds to exactly 1. So the computed $A^\top A$ is $\begin{pmatrix} 1 & 1 \\ 1 & 1 \end{pmatrix}$ , singular, and solve stops with an error.

lstsq still finds $\hat{\mathbf{x}} = (1, 1)$ for $\mathbf{b} = (2, 10^{-8}, 10^{-8})$ .

Rule: In code

\texttt{np.linalg.lstsq(A, b, rcond=None)}

Use it for least squares rather than solving $A^\top A\,\hat{\mathbf{x}} = A^\top\mathbf{b}$ , which squares the condition number. rcond=None sets how small a singular value counts as zero; it is the setting NumPy recommends.

A^\top A

import numpy as np years = np.arange(2015, 2025, dtype=float) A = np.column_stack([np.ones(10), years]) centred = np.column_stack([np.ones(10), years - years.mean()]) print(f"{np.linalg.cond(A):.1e}", f"{np.linalg.cond(A.T @ A):.1e}", round(np.linalg.cond(centred), 4)) cos = np.ones(10) @ years / (np.sqrt(10) * np.linalg.norm(years)) print(round(np.degrees(np.arccos(cos)), 2)) eps = 1e-8 L = np.array([[1, 1], [eps, 0], [0, eps]]) b = np.array([2, eps, eps]) print(L.T @ L) try: np.linalg.solve(L.T @ L, L.T @ b) except np.linalg.LinAlgError as err: print("LinAlgError:", err) print(np.round(np.linalg.lstsq(L, b, rcond=None)[0], 4))

⌘/Ctrl + Enter runs the cell

When the columns are dependent

If the columns of $A$ are dependent (lesson 3.5.1), $A^\top A$ has no inverse. The projection $\mathbf{p}$ is still unique: it is the nearest point of the column space. But many $\hat{\mathbf{x}}$ give that same $\mathbf{p}$ .

Fit the five points with a third column that is twice the second, $y = c + m_1x + m_2(2x)$ . Every choice with $c = 1$ and $m_1 + 2m_2 = 1.2$ gives the fitted values $(1, 2.2, 3.4, 4.6, 5.8)$ : $(1, 1.2, 0)$ , $(1, 0, 0.6)$ and $(1, 0.24, 0.48)$ among them.

lstsq does not fail here. It reports rank 2 instead of 3, returns the shortest $\hat{\mathbf{x}}$ , here $(1, 0.24, 0.48)$ , and leaves its sum of squares empty. Read the rank before reading the answer.

The same happens when every $x_i$ is the same: the column of $x$ values is a multiple of the column of ones. For points $(2, 1)$ , $(2, 3)$ , $(2, 5)$ , a best line would be vertical, which has no slope (lesson 1.2.1). lstsq reports rank 1 and returns $(0.6, 1.2)$ , which means nothing as an intercept and slope.

Rule: Dependent columns

\begin{gathered} \mathbf{p} \text{ is unique;} \\ \hat{\mathbf{x}} \text{ is not} \end{gathered}

When the rank is below the number of columns, many $\hat{\mathbf{x}}$ give $\mathbf{p}$ . Each is one of them plus a vector of the null space (lesson 3.5.1). Check the rank that lstsq returns.

\mathbf{1}

import numpy as np x = np.array([0, 1, 2, 3, 4], dtype=float) y = np.array([1, 2, 4, 4, 6], dtype=float) D = np.column_stack([np.ones(5), x, 2 * x]) x_hat, ss, rank, _ = np.linalg.lstsq(D, y, rcond=None) print(np.round(x_hat, 4), ss, rank) for guess in [(1, 1.2, 0), (1, 0, 0.6), tuple(x_hat)]: print(np.round(guess, 4), np.round(D @ np.array(guess), 4)) same = np.column_stack([np.ones(3), [2.0, 2.0, 2.0]]) sol, _, rank, _ = np.linalg.lstsq(same, np.array([1.0, 3.0, 5.0]), rcond=None) print(np.round(sol, 4), rank)

⌘/Ctrl + Enter runs the cell

Case study 3.5.5 uses least squares twice: for the minimum-variance hedge, and for the closest copy of a payoff that cannot be copied exactly. Module 5.10 adds the statistics: how far $\hat{\mathbf{x}}$ may be from the truth.

Common mistakes

Least squares counts squared misses, not equations met. $(2, 1.5)$ meets two of the three equations exactly, but its total, 0.5, is larger than the least-squares 0.2857.
Residuals are vertical gaps, misses in $y$ . They are not the distances from the points to the line at right angles.
Keep the column of ones. Without it the line must pass through the origin: $y = mx$ gives $m = \frac{46}{30} \approx 1.5333$ and a sum of squares of 2.4667, not 0.8.
The order of the columns is the order of $\hat{\mathbf{x}}$ . With the ones first, the intercept comes first; np.polyfit lists the slope first.
Do not form $A^\top A$ for real data. It squares the condition number; use lstsq, and centre data such as calendar years.
Check the rank. With dependent columns lstsq still returns an answer, but not the only one.

Least squares

0 of 11 correct

1.
For $\mathbf{x} = (2, 1)$ in the system above, $T\mathbf{x} - \mathbf{b} = (0, 0, -1)$ . What is $\lVert T\mathbf{x} - \mathbf{b} \rVert^2$ ?
2.
For $T = \begin{pmatrix} 1 & 1 \\ 1 & -1 \\ 1 & 2 \end{pmatrix}$ , find the entry in row 2, column 2 of $T^\top T$
3.
Fit $y = c + mx$ to $(0, 1)$ , $(1, 1)$ , $(2, 4)$ . Find $m$
4.
For the same three points, find $c$
5.
Fit $y = c + mx$ to $(-1, 1)$ , $(0, 2)$ , $(1, 6)$ . Find $m$
6.
For that fit, $y = 3 + 2.5x$ , what is the sum of squared residuals?
7.
The best constant $c$ for the data 2, 3, 7 solves $A\mathbf{x} = \mathbf{y}$ with $A$ a column of ones. What is $c$ ?
8.
Fit $y = mx$ , a line through the origin, to $(1, 2)$ and $(2, 3)$ . Find $m$
9.
A line fit with an intercept has five residuals. Four of them are 0.3, $-0.5$ , 0.1 and 0.4. What is the fifth?
10.
$\operatorname{cond}(A) = 1000$ . About how large is $\operatorname{cond}(A^\top A)$ ? Give a power of 10
11.
$A$ is $5 \times 3$ and its third column is the sum of the first two. What rank does lstsq report?

How to type answers

Powers with ^: x^7, 2^(1/3), a^-2.
Multiply with * or by writing side by side: 3x^2, 2(x+1).
ln(x) is the natural log, log(x, 2) is log base 2, e is Euler's number, sqrt(x) is a square root.
sin(x), cos(x) and tan(x) take radians; write pi for π and sin(x)^2 for sin²x.
Your first answer to each question is saved: it counts towards your score and brings the question back for review.
Exact answers like ln(20)/ln(3) are best. Decimals need about 4 significant figures.

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

Least squares

A system with no solution

The best answer

Rule: Least squares

The normal equations

Rule: The normal equations

Fitting a line to points

The slope as a projection

Rule: Least-squares line

A noisy line

Solve with lstsq, not with AᵀA

Rule: In code

When the columns are dependent

Rule: Dependent columns

Common mistakes

Least squares

Summary card

Free further reading