Orthogonality, projection and least squares

Suppose you record the number of requests hitting a server and its response time over a day, and plot them as points. The points roughly form a band rising to the upper right, but they do not sit exactly on one line, because measurements always carry noise. Still, you want to know "when requests go up by 100, roughly how many ms does the response time go up?" You need the one line that passes through the points as well as possible.

Through the eyes of Lesson 8, this is a strange system of equations. Each point is one equation, so 100 points make 100 equations, yet there are only two unknowns, the slope and the intercept. A system like this almost always has no solution. In this lesson you will learn how to find "the closest answer" when there is no solution: the method of least squares. Linear regression in statistics, the most basic model in machine learning, sensor calibration and trend lines on charts are all this method.

The key is the shadow (projection) from Lesson 3 and the idea of "perpendicular". The closest point is always at the foot of the perpendicular.

13-1. Orthogonal: the dot product is 0 (5 min)

In Lesson 3 you learned that two vectors are perpendicular when their dot product is 0. In linear algebra, perpendicular is called orthogonal. For example, $(2, 1) \cdot (-1, 2) = -2 + 2 = 0$, so the two vectors are orthogonal.

A set of vectors that are orthogonal to each other and all have length 1 is called orthonormal. The standard basis î, ĵ is the classic example, and a rotated copy of it is orthonormal too.

An orthonormal basis is convenient because finding coordinates takes just one dot product. In Lesson 10 you had to solve $B\mathbf{c} = \mathbf{v}$ to get the new coordinates. But if $\mathbf{q}_1, \mathbf{q}_2$ are orthonormal, the coordinates of $\mathbf{v}$ are simply $c_1 = \mathbf{q}_1 \cdot \mathbf{v}$ and $c_2 = \mathbf{q}_2 \cdot \mathbf{v}$. The length of the shadow on each axis is the coordinate. We will put this convenience straight to use in 13-5.

13-2. Projecting onto a line (6 min)

In Lesson 3 you learned that the length of the shadow $\mathbf{w}$ casts on $\mathbf{u}$ is $\dfrac{\mathbf{u} \cdot \mathbf{w}}{\lVert \mathbf{u} \rVert}$. Now let's find the shadow as a vector, not a length. The shadow points along $\mathbf{u}$, so it is some multiple of $\mathbf{u}$.

This is called the projection of $\mathbf{w}$ onto $\mathbf{u}$. The denominator is $\mathbf{u} \cdot \mathbf{u} = \lVert\mathbf{u}\rVert^2$ rather than $\lVert\mathbf{u}\rVert$ because the shadow length was multiplied by the unit vector $\mathbf{u}/\lVert\mathbf{u}\rVert$ (dividing by $\lVert\mathbf{u}\rVert$ once more).

Example. If $\mathbf{u} = (2, 1)$ and $\mathbf{w} = (1, 3)$, then $\mathbf{u} \cdot \mathbf{w} = 2 + 3 = 5$ and $\mathbf{u} \cdot \mathbf{u} = 4 + 1 = 5$, so $\mathbf{p} = 1 \cdot (2, 1) = (2, 1)$.

Let's call the leftover part $\mathbf{e} = \mathbf{w} - \mathbf{p} = (-1, 2)$ the error. $\mathbf{e} \cdot \mathbf{u} = -2 + 2 = 0$: the error is orthogonal to the line. This is the key property of projection.

A picture shows why right away. The shortest way from a point to a line is straight down at a right angle. Going at a slant means riding the hypotenuse of a right triangle, which is always longer (Pythagoras). Least squares is this idea, "perpendicular means closest", extended from a line to a larger space.

13-3. An Ax = b that cannot be solved: fitting a line to points (7 min)

Let's start with a small example. You want to fit a line $y = mx + c$ to the four points $(0, 1)$, $(1, 3)$, $(2, 2)$, $(3, 4)$. Each point gives one equation saying "the line must pass through this point".

Column 1 of the matrix $X$ holds the points' $x$ values, and column 2 is all 1s. The 1s in column 2 are the slots that get multiplied by the intercept $c$. There are four equations and two unknowns. The first two equations give $c = 1$ and $m = 2$, but plugging those into the third gives $2 \cdot 2 + 1 = 5 \neq 2$. There is no solution.

The column picture from Lesson 9 makes the reason clear. The places $X\mathbf{c}$ can reach (the column space) form a single plane spanned by the two columns of $X$. $\mathbf{y}$ is a point in 4-dimensional space that does not lie on that plane. So you cannot land on it exactly, and the best you can do is go to the point on the plane closest to $\mathbf{y}$. By 13-2, that point is the projection of $\mathbf{y}$ onto the plane.

The vertical gap from each point to the line is called a residual.

$\lVert \mathbf{r} \rVert^2$ is called the sum of squared residuals, and choosing the $\hat{\mathbf{c}}$ that makes it as small as possible is least squares. Making the distance as short as possible = making the squared distance as small as possible = making the sum of squared residuals as small as possible: they all say the same thing.

There are two reasons to square the residuals instead of just adding them. If you just add them, the points above the line (+) and the points below (−) cancel each other out. And squaring matches the Pythagorean "distance", so all the projection and orthogonality tools you have learned so far still apply.

The lab below has 6 black data points. The thick blue (indigo) line is the least-squares line, and the crimson dashed lines are the residuals (vertical gaps) at each point. The line's equation and the sum of squared residuals appear on the right. At first it shows "y = 0.717x + 0.393" with a sum of squared residuals of 1.007.

Turn on 'My line' and a teal dashed line appears along with its own sum of squared residuals. However hard you try, the message "The least-squares line wins. No line can make this sum smaller" does not change, until your line almost overlaps the blue one and "Almost equal! You found the least-squares line" appears. The least-squares line is, as the name says, the line with the smallest sum of squared residuals, so you cannot beat it, only tie.

Also notice how much the line bends when you drag one point far away. Because the residuals are squared, a single point far off inflates the whole sum, and the line gets dragged toward it to appease it. Remember that least squares is sensitive to outliers.

Finally, a residual is the vertical gap, not the perpendicular distance to the line. The problem is "predict y when you know x", so only the error in the y direction is measured.

13-4. The normal equations: the residual is orthogonal to every column (7 min)

Extend the property from 13-2 to a plane and you can compute the answer. If you went from $\mathbf{y}$ to the closest point on the plane (the column space), the leftover residual $\mathbf{r}$ is orthogonal to the plane. Being orthogonal to the plane means being orthogonal to every column that spans it. Writing the dot products with all the columns at once gives $X^\mathsf{T}\mathbf{r}$, so we get this.

These are called the normal equations. "Normal" here does not mean "ordinary"; it means "perpendicular". $X$ is 4×2, so it is not square, but $X^\mathsf{T}X$ is a 2×2 square matrix, so you can solve it as in Lesson 8.

Let's solve it by hand. For the four points from 13-3, the calculation goes like this.

Solving $14m + 6c = 19$ and $6m + 4c = 10$ gives $m = 0.8$ and $c = 1.3$. The best-fitting line is $y = 0.8x + 1.3$.

Check. The line's values at $x = 0, 1, 2, 3$ are $1.3, 2.1, 2.9, 3.7$, and the residuals are $\mathbf{r} = (-0.3,\ 0.9,\ -0.9,\ 0.3)$.

• Dot product with column 2 (all 1s): $-0.3 + 0.9 - 0.9 + 0.3 = 0$. The residuals add up to 0, so the amounts above and below the line balance out.
• Dot product with column 1 (the $x$ values): $0 + 0.9 - 1.8 + 0.9 = 0$.
• Sum of squared residuals: $0.09 + 0.81 + 0.81 + 0.09 = 1.8$.

The residual is orthogonal to both columns, so this line really is the closest answer.

13-5. QR decomposition: switching to an orthogonal basis makes it easy (6 min)

The normal equations are neat, but they have a weakness when solved on a computer. The moment you form $X^\mathsf{T}X$, numerical errors tend to grow, and the more the columns point in similar directions, the worse it gets. The property from 13-1, "in an orthonormal basis, coordinates take one dot product", lets you avoid this problem.

The Gram–Schmidt process takes any columns and turns them into an orthonormal basis that spans the same plane.

1. Scale the first column $\mathbf{a}_1$ to length 1: $\mathbf{q}_1 = \mathbf{a}_1 / \lVert\mathbf{a}_1\rVert$.
2. Subtract the shadow in the $\mathbf{q}_1$ direction from the second column $\mathbf{a}_2$: $\mathbf{b} = \mathbf{a}_2 - (\mathbf{q}_1 \cdot \mathbf{a}_2)\,\mathbf{q}_1$. Like the error in 13-2, the leftover $\mathbf{b}$ is orthogonal to $\mathbf{q}_1$.
3. Scale $\mathbf{b}$ to length 1: $\mathbf{q}_2 = \mathbf{b} / \lVert\mathbf{b}\rVert$.

If there are more columns, you keep going: subtract the shadows along all the $\mathbf{q}$s made so far, then scale to length 1. Call the matrix with these $\mathbf{q}$s as its columns $Q$, and you can write the original matrix like this.

The columns of $Q$ are orthonormal, and $R$ is an upper triangular matrix, with zeros below the diagonal. $R$ records "how the original columns were combined from the $\mathbf{q}$s"; $\mathbf{a}_1$ is made from $\mathbf{q}_1$ alone and $\mathbf{a}_2$ from $\mathbf{q}_1, \mathbf{q}_2$, which gives the upper triangular shape. This is called the QR decomposition.

Put $X = QR$ into the normal equations, and $Q^\mathsf{T}Q = I$ shrinks the equation.

$Q^\mathsf{T}\mathbf{y}$ is "the orthonormal coordinates of $\mathbf{y}$", just a few dot products. $R$ is upper triangular, so you solve one row at a time starting from the bottom and you are done (the bottom equation has only one unknown). You never form $X^\mathsf{T}X$, so errors grow less. You will see how much less, in numbers, in 13-6.

13-6. Least squares with NumPy (12 min)

Reading the lab's NumPy panel

On the lab's starting screen, the NumPy panel shows this code.

>>> import numpy as np
>>> x = np.array([-3.5, -2, -1, 0.5, 2, 3.5])
>>> y = np.array([-2, -1.5, 0, 0.5, 2.5, 2.5])
>>> X = np.column_stack([x, np.ones_like(x)])
>>> coef, res, rank, sv = np.linalg.lstsq(X, y)
>>> np.round(coef, 6)
array([0.716934, 0.393078])
>>> np.round(res, 6)
array([1.007417])
>>> int(rank)
2

Let's go line by line. x and y are the x and y coordinates of the 6 points. Some are written as integers like -2, but since there is at least one decimal like 3.5, the whole array is float64 (Lesson 1).

np.ones_like(x) means "an array of all 1s with the same shape and dtype as x". It makes one of the intercept 1s from 13-3 for each point. Stand the two arrays up as columns with np.column_stack from Lesson 4, and you get the matrix X from 13-3.

>>> x.dtype
dtype('float64')
>>> np.ones_like(x)
array([1., 1., 1., 1., 1., 1.])
>>> X
array([[-3.5,  1. ],
       [-2. ,  1. ],
       [-1. ,  1. ],
       [ 0.5,  1. ],
       [ 2. ,  1. ],
       [ 3.5,  1. ]])
>>> X.shape
(6, 2)

X has 6 rows and 2 columns: 6 equations, 2 unknowns.

np.linalg.lstsq(X, y) finds the least-squares solution of Xc ≈ y. The name is short for least squares. There are four results, so you receive them under four names.

Let's also look at the unrounded values and receive the slope and intercept separately.

>>> coef
array([0.71693449, 0.39307787])
>>> m, c = coef
>>> res
array([1.00741656])
>>> sv
array([5.81025009, 2.44765069])

There is a reason the panel shows rank as int(rank). Printed as is, rank comes out as a different kind of integer depending on the operating system: np.int32(2) on Windows, np.int64(2) on Linux. The value is 2 either way, so converting it to a Python integer with int(…) makes it look the same everywhere.

Checking that the residual is orthogonal to the columns

Let's check the property from 13-4 directly. X @ coef gives the values on the line (the projected point), and y - X @ coef is the residual.

>>> r = y - X @ coef
>>> r
array([ 0.11619283, -0.4592089 ,  0.32385661, -0.25154512,  0.67305315,
       -0.40234858])
>>> r @ r
np.float64(1.007416563658838)
>>> np.allclose(X.T @ r, 0)
True

r @ r is the sum of the squared components, so it is the sum of squared residuals, the same as res. X.T @ r computes the dot products of the residual with both columns at once, and it comes out as 0. It may not be an exact 0, with tiny residue left over, so we compared with np.allclose.

Let's compare with "my line" too. In the lab, 'My line' starts as y = x/6 + 0.5, through the points (−3, 0) and (3, 1). Write a small function for the sum of squared residuals and compare.

>>> def sse(m, c):
...     r = y - (m * x + c)
...     return r @ r
>>> sse(m, c)
np.float64(1.007416563658838)
>>> sse(1/6, 0.5)
np.float64(11.354166666666668)
>>> sse(0.7, 0.4)
np.float64(1.0175)

Even (0.7, 0.4), a pretty good fit by eye, gives 1.0175, a little more than the least-squares line's 1.0074.

Solving the normal equations with solve

Solve XᵀX ĉ = Xᵀy from 13-4 with np.linalg.solve from Lesson 8 and see whether it matches lstsq.

>>> X.T @ X
array([[33.75, -0.5 ],
       [-0.5 ,  6.  ]])
>>> X.T @ y
array([24.,  2.])
>>> np.linalg.solve(X.T @ X, X.T @ y)
array([0.71693449, 0.39307787])
>>> np.allclose(np.linalg.solve(X.T @ X, X.T @ y), coef)
True

Solve the four points from 13-3 the same way to confirm the y = 0.8x + 1.3 you found by hand.

>>> xs = np.array([0, 1, 2, 3])
>>> ys = np.array([1, 3, 2, 4])
>>> Xs = np.column_stack([xs, np.ones_like(xs)])
>>> Xs.T @ Xs
array([[14,  6],
       [ 6,  4]])
>>> Xs.T @ ys
array([19, 10])
>>> np.linalg.solve(Xs.T @ Xs, Xs.T @ ys)
array([0.8, 1.3])
>>> np.linalg.lstsq(Xs, ys)[1]
array([1.8])

lstsq(…)[1] takes out only the second of the four results (res). It matches the sum of squared residuals 1.8 you found by hand.

QR decomposition: np.linalg.qr

>>> Q, R = np.linalg.qr(X)
>>> Q.shape, R.shape
((6, 2), (2, 2))
>>> np.round(Q, 4)
array([[-0.6025, -0.3873],
       [-0.3443, -0.3964],
       [-0.1721, -0.4024],
       [ 0.0861, -0.4115],
       [ 0.3443, -0.4206],
       [ 0.6025, -0.4297]])
>>> np.round(R, 4)
array([[ 5.8095, -0.0861],
       [ 0.    , -2.448 ]])
>>> np.allclose(Q.T @ Q, np.eye(2))
True
>>> np.allclose(Q @ R, X)
True
>>> np.linalg.solve(R, Q.T @ y)
array([0.71693449, 0.39307787])

Q has the same 6 rows and 2 columns as X, and its two columns are orthonormal, so Q.T @ Q is the 2×2 identity matrix. R is upper triangular, with zeros below the diagonal. Solving R ĉ = Qᵀy gives the same answer.

The many negative numbers in Q and R come from a freedom in the signs. Flip the sign of a column of Q together with the same row of R and the product QR stays the same, so NumPy picks whichever is convenient to compute (the same story as the signs of eigenvectors in Lesson 11). Follow the Gram–Schmidt process from 13-5 yourself and compare with Q: they are the same except for the signs.

>>> a1, a2 = X[:, 0], X[:, 1]
>>> q1 = a1 / np.linalg.norm(a1)
>>> b = a2 - (q1 @ a2) * q1
>>> q2 = b / np.linalg.norm(b)
>>> round(float(q1 @ q2), 10)
0.0
>>> np.allclose(np.abs(Q), np.abs(np.column_stack([q1, q2])))
True

np.abs takes the absolute value of each element, erasing the signs.

There is also a way to see in numbers what 13-5 said, that "forming XᵀX tends to make errors grow". np.linalg.cond gives a value called the condition number, which says "by at most how many times a small error in the input can grow in the answer".

>>> round(float(np.linalg.cond(X)), 3)
2.374
>>> round(float(np.linalg.cond(X.T @ X)), 3)
5.635

2.374² ≈ 5.635: the condition number of XᵀX is the square of the condition number of X. In this example both are small, so there is no problem, but if X's condition number were 10⁸, XᵀX's would be 10¹⁶, eating up all the precision of floating point (about 16 digits). That is why lstsq does not solve the normal equations directly and uses a more stable method (based on SVD) instead. The condition number is a central topic of the next course, numerical analysis.

np.polyfit: a trend line in one line

np.polyfit(x, y, degree), which fits a polynomial to data, solves the same least-squares problem under the hood. Degree 1 is a straight line, and the result is ordered from the highest degree down: [slope, intercept].

>>> np.polyfit(x, y, 1)
array([0.71693449, 0.39307787])
>>> np.allclose(np.polyfit(x, y, 1), coef)
True

Curves work the same way. To fit y = ax² + bx + c, just add one more column of x² to X. In the unknowns a, b and c it is still a "linear" equation, so linear algebra solves it.

>>> xq = np.array([-2, -1, 0, 1, 2])
>>> yq = np.array([4.1, 1.2, 0.1, 0.9, 4.2])
>>> np.linalg.lstsq(np.column_stack([xq**2, xq, np.ones_like(xq)]), yq)[0]
array([ 1.02142857, -0.01      ,  0.05714286])
>>> np.polyfit(xq, yq, 2)
array([ 1.02142857, -0.01      ,  0.05714286])

y ≈ 1.02x² − 0.01x + 0.06: we found a parabola that is almost y = x².

13-7. Common mistakes (3 min)

Passing a non-square X to solve. solve only handles square systems, with as many equations as unknowns (Lesson 8). When there are more equations, you have to use lstsq.

>>> np.linalg.solve(X, y)
Traceback (most recent call last):
  ...
numpy.linalg.LinAlgError: Last 2 dimensions of the array must be square

Leaving out the column of 1s. Pass only x, without the np.ones_like(x) column, and you only find y = mx, a line through the origin. You have pinned the intercept to 0, so the slope changes too.

>>> np.linalg.lstsq(x.reshape(-1, 1), y)[0]
array([0.71111111])

x.reshape(-1, 1) stands x up as a matrix with 6 rows and 1 column (-1 means "work out the count yourself"). The slope came out as 0.711, not 0.717.

When res is empty. If the rank of X is less than the number of columns (for example, every point has the same x), res comes back as an empty array. In the lab, this is the case where you gather all 6 points on one vertical line and the warning "All x values are the same, so no single line is defined" appears. Even then lstsq returns an answer without any error (the shortest of the possible answers), but that does not mean a single line has been determined. If res is empty, check rank first.

>>> xv = np.array([1.0, 1, 1, 1, 1, 1])
>>> coef2, res2, rank2, sv2 = np.linalg.lstsq(np.column_stack([xv, np.ones_like(xv)]), y)
>>> res2
array([], dtype=float64)
>>> int(rank2)
1

Exercises (4 min)

1. Find the projection $\mathbf{p}$ of $\mathbf{w} = (3, 1)$ onto $\mathbf{u} = (1, 1)$ and the error $\mathbf{w} - \mathbf{p}$, and check that the error is orthogonal to $\mathbf{u}$.
2. Find the least-squares line for the points $(1, 2)$, $(2, 3)$, $(3, 5)$, $(4, 6)$ in two ways, with lstsq and with polyfit, and find the sum of squared residuals too.
3. Use X.T @ r to check that the residual from question 2 is orthogonal to both columns of $X$. What do the residuals add up to?

Solutions

Question 1: $\mathbf{u} \cdot \mathbf{w} = 4$ and $\mathbf{u} \cdot \mathbf{u} = 2$, so $\mathbf{p} = 2\,(1, 1) = (2, 2)$, the error is $(1, -1)$, and $(1, -1) \cdot (1, 1) = 0$.

Question 2: $y = 1.4x + 0.5$, the residuals are $(0.1, -0.3, 0.3, -0.1)$, and the sum of squared residuals is $0.01 + 0.09 + 0.09 + 0.01 = 0.2$.

Question 3: column 2 is all 1s, so "the dot product with column 2 = the sum of the residuals", which is 0.

>>> u = np.array([1, 1])
>>> w = np.array([3, 1])
>>> p = (u @ w) / (u @ u) * u
>>> p
array([2., 2.])
>>> (w - p) @ u
np.float64(0.0)
>>> xe = np.array([1, 2, 3, 4])
>>> ye = np.array([2, 3, 5, 6])
>>> Xe = np.column_stack([xe, np.ones_like(xe)])
>>> coef, res, rank, sv = np.linalg.lstsq(Xe, ye)
>>> coef
array([1.4, 0.5])
>>> res
array([0.2])
>>> np.polyfit(xe, ye, 1)
array([1.4, 0.5])
>>> re = ye - Xe @ coef
>>> re
array([ 0.1, -0.3,  0.3, -0.1])
>>> np.allclose(Xe.T @ re, 0)
True

The line that computes p is the formula from 13-2 written out directly. (u @ w) / (u @ u) is a single number (the factor), and it multiplies the vector u.

In the next lesson you will learn the singular value decomposition (SVD), which splits every matrix into "rotate, stretch, rotate again". It is what the sv returned by today's lstsq really is, and it is also the stable method lstsq uses under the hood.