Projection, Orthogonality and Least Squares

Context

Projecting one vector onto another finds the point on the second vector's line that sits closest to the first, and the piece left over, the residual, meets that line at a right angle. Least squares makes the same move with a whole table of features: it picks the weights whose predictions land closest to the real targets.

Every linear regression (a model that predicts a number as a weighted sum of the features plus a bias) computes this. Fitting price on area across a set of listings, giving the model a column of ones so it learns an intercept, 2.3.1's bias in w . x + b, and reading the residuals a fit leaves behind are all projection at work. The textbook formula for a regression's weights is the same computation written as a system of equations, and tutorial code types it out in full.

The failures look like working code. A fit with no column of ones is forced through zero price at zero area, and nothing warns. A residual tested with == 0 fails on float noise for a correct fit. The second value np.linalg.lstsq returns gets read as the residuals when it is a single number. Code that inverts a matrix (computes the matrix that undoes a multiplication by it) for the weights matches on a tidy example and misses by far more on a real column such as a construction year. The previous lesson measured the distance between two points, this lesson finds the point at the smallest distance.

Goal

• Compute the projection of one vector onto another, and its scalar projection, with @ and np.linalg.norm.
• Verify that a residual is orthogonal to what it was projected onto, with a test that survives float noise and large numbers.
• Fit price on area with an intercept column through np.linalg.lstsq, unpacking its four return values correctly.
• Solve the normal equations with np.linalg.solve and confirm they reproduce the weights lstsq found.
• Explain precisely why np.linalg.inv is the wrong tool for solving a system of equations.

Prerequisites

None of this is required, the lesson stands on its own. Both are here if you want the geometry drawn on a board, or the function's reference open, before you start.

• YouTube - "15. Projections onto Subspaces" by MIT OpenCourseWare, a lecture by Gilbert Strang. It projects a vector onto a line, then onto a plane, and arrives at the normal equations, the path this lesson takes in code.
• Official docs - numpy.linalg.lstsq, the NumPy 2.5 reference page. The four return values, and when the residuals array comes back empty.

Table of Contents

• Projecting One Vector onto Another
• The Residual Is Orthogonal
• Least Squares with np.linalg.lstsq
• The Normal Equations

Theory

Projecting One Vector onto Another

What. The projection of a onto b is the vector on the line through b that lies closest to a, in plain text proj = (a . b / b . b) b. The scalar projection, a . b / |b|, is that vector's signed length. Why. A model with one feature and no intercept can only predict multiples of the feature's column, so its best fit is a projection: the multiple of the area column that lands closest to the price column. When. Whenever one vector has to be approximated by a multiple of another: a one-feature fit through the origin, how far a vector reaches along a direction, or the part of a target that one feature accounts for.

The fraction a . b / b . b counts how many copies of b to take, and multiplying b by it lands on b's line. Two vectors small enough to draw keep every number checkable by eye.

import numpy as np

a = np.array([1.0, 3.0])
b = np.array([4.0, 2.0])

# a . b / b . b, here 10 / 20
coef = (a @ b) / (b @ b)
proj = coef * b  # that many copies of b
print(coef)      # 0.5
print(proj)      # [2. 1.]

The scalar projection measures the same reach as a length along b. It uses 2.3.3's np.linalg.norm for |b|, and its sign is the sign of the dot product from 2.3.1, positive when a points along b and negative when it points away.

import numpy as np

a = np.array([1.0, 3.0])
b = np.array([4.0, 2.0])

# a . b / |b|, a signed length measured along b
scalar = (a @ b) / np.linalg.norm(b)
print(np.round(scalar, 4))  # 2.2361
# 2.2361, the length of proj = [2. 1.] itself
print(np.round(np.linalg.norm(0.5 * b), 4))
# -2.2361, a turned around points away from b
print(np.round(((-a) @ b) / np.linalg.norm(b), 4))

Nothing in the formula needs two entries. A price column with one entry per listing is a vector, and so is the area column beside it, and the c below comes out as a price per square meter for a line forced through the origin.

import numpy as np

# five listings: area in square meters, price in thousands
area = np.array([50.0, 60.0, 80.0, 100.0, 120.0])
price = np.array([198.0, 210.0, 285.0, 332.0, 405.0])

# the multiple of area that lands closest to price
c = (price @ area) / (area @ area)
fitted = c * area           # the model price = c * area
print(np.round(c, 4))       # 3.4444
print(np.round(fitted, 2))  # [172.22 206.67 275.56 344.44 413.33]

Closest here is 2.3.3's Euclidean distance between the price vector and the vector of predictions, measured in a space with one axis per listing. That space cannot be drawn, and @ does not need it drawn.

DON'T:

• DON'T put np.linalg.norm(b) in the denominator when you want the projection vector. (a @ b) / np.linalg.norm(b) * b returns [8.94427191 4.47213595], |b| times longer than [2. 1.], because the scalar projection is a length, not a count of copies of b.
• DON'T project onto a column of zeros. b @ b is 0.0, and the division fails the way 2.3.1's cosine of an all-zero vector did, nan with the same RuntimeWarning, so a feature that was never filled in turns every number computed from it into nan.

The Residual Is Orthogonal

What. The residual, r = a - proj, is the part of a that b's line cannot reach, and r . b is zero: the residual is orthogonal to b, 2.3.1's word for perpendicular. Why. The right angle is what makes the projection the closest point. Every other point on the line sits farther from a, so the residual's length is the shortest distance from a to the line. When. Every time you read a residual, the part of the target the features could not explain, or confirm that a fit is the best its features allow.

Subtracting the projection from a leaves the residual, and for these small vectors its dot product with b is exactly zero.

import numpy as np

a = np.array([1.0, 3.0])
b = np.array([4.0, 2.0])
proj = (a @ b) / (b @ b) * b

r = a - proj  # what b's line cannot reach
print(r)      # [-1.  2.]
print(r @ b)  # 0.0

That the right angle marks the closest point can be measured. Every point on the line is a multiple t * b, and the block below measures the distance from a to five of them, with a broadcast from 2.2.4 and a row-wise norm from 2.3.3.

import numpy as np

a = np.array([1.0, 3.0])
b = np.array([4.0, 2.0])

ts = np.array([0.0, 0.25, 0.5, 0.75, 1.0])  # steps along b
points = ts[:, np.newaxis] * b              # (5, 1) * (2,) is (5, 2)
dist = np.linalg.norm(a - points, axis=1)   # a's distance to each row
print(np.round(dist, 4))                    # [3.1623 2.5    2.2361 2.5    3.1623]
print(ts[dist.argmin()])                    # 0.5
print(points)

[[0.  0. ]
 [1.  0.5]
 [2.  1. ]
 [3.  1.5]
 [4.  2. ]]

The shortest distance, 2.2361, belongs to t = 0.5, the coef the projection computed, and it is the residual's own length. The points on either side are farther away by the same amount, 2.5 a quarter step in each direction.

On listings the zero stops being exact. The price per square meter c is 31/9, which no float holds exactly, so the residual carries rounding and its dot product with area misses zero by a sliver.

import numpy as np

area = np.array([50.0, 60.0, 80.0, 100.0, 120.0])
price = np.array([198.0, 210.0, 285.0, 332.0, 405.0])
c = (price @ area) / (area @ area)

r = price - c * area            # one residual per listing
print(np.round(r, 2))           # [ 25.78   3.33   9.44 -12.44  -8.33]
print(r @ area)                 # -5.9117155615240335e-12
print(r @ area == 0)            # False
print(np.isclose(r @ area, 0))  # True
print(np.round(r @ r, 2))       # 989.11

The sliver's digits depend on the machine, so yours can differ, and np.isclose, from 2.3.1, accepts it where == does not. The last line squares the residual's length: r @ r adds up each listing's squared residual, so the closest point is also the price per square meter with the smallest sum of squared residuals, the quantity least squares is named after.

Against zero, np.isclose passes anything within 1e-08, a fixed allowance that ignores how large the numbers behind the dot product were. Price the same listings in plain currency units and measure their area in square feet, 10.7639 to the square meter, and the rounding outgrows that allowance, while 2.3.1's cosine, r . b / (|r| |b|), carries no units and still sits at zero.

import numpy as np

# the same listings, square feet and plain currency units
sqft = np.array([50.0, 60.0, 80.0, 100.0, 120.0]) * 10.7639
price = np.array([198.0, 210.0, 285.0, 332.0, 405.0]) * 1000

r = price - (price @ sqft) / (sqft @ sqft) * sqft
print(np.isclose(r @ sqft, 0))  # False
cos = (r @ sqft) / (np.linalg.norm(r) * np.linalg.norm(sqft))
print(np.isclose(cos, 0))       # True

The right angle lives in the space with one axis per listing. Plotted as points of area against price, the same five residuals are vertical gaps between each point and the line, gaps that meet the line at a slant.

DON'T:

• DON'T test a residual with == 0. A correct projection printed False above, because the float holding c is rounded, and np.isclose or the cosine is the test.
• DON'T look for the right angle in a scatter plot. Area against price shows each residual as a vertical gap that meets the line at a slant, and the right angle exists only between whole vectors, one axis per listing.

Least Squares with np.linalg.lstsq

What. Least squares chooses the weights w whose predictions X @ w land closest to the targets y, which is the same as the smallest sum of squared residuals. np.linalg.lstsq(X, y) computes them. Why. With more listings than weights, X w = y has an exact solution only when every listing sits exactly on the model, and five listings that do not lie on one line cannot all sit on one. Least squares gives the closest fit the columns allow, the fit every linear regression computes. When. Every time a linear model is fitted on a design matrix: price on area with an intercept, or price on several features at once.

The model price = w0 + w1 * area needs a column for w0 to multiply, and a column of ones is that column, so each row's prediction is w0 * 1 + w1 * area. Lesson 2.2.1's np.column_stack builds the design matrix from the ones and the areas.

import numpy as np

area = np.array([50.0, 60.0, 80.0, 100.0, 120.0])
price = np.array([198.0, 210.0, 285.0, 332.0, 405.0])

# a column of ones beside the column of areas
X = np.column_stack([np.ones(len(area)), area])
print(X.shape)  # (5, 2)
print(X)

[[  1.  50.]
 [  1.  60.]
 [  1.  80.]
 [  1. 100.]
 [  1. 120.]]

np.linalg.lstsq returns a tuple of four, unpacked here with Module 1's tuple unpacking: the weights, the sum of squared residuals inside a one-element array, the rank of X (the number of independent columns, Lesson 2.3.5's subject), and the singular values of X (numbers describing how strongly X stretches space in each direction, which 2.3.5's condition number is built from).

import numpy as np

area = np.array([50.0, 60.0, 80.0, 100.0, 120.0])
price = np.array([198.0, 210.0, 285.0, 332.0, 405.0])
X = np.column_stack([np.ones(len(area)), area])

w, ssr, rank, sv = np.linalg.lstsq(X, price)
print(np.round(w, 4))   # [40.  3.]
print(ssr)              # [278.]
print(rank)             # 2
print(np.round(sv, 4))  # [192.1056   0.6666]

Code written for older NumPy passes rcond=None to this call. In NumPy 2.5 that is already the default, so the line above makes the same call.

The predictions X @ w are the projection of the price vector onto everything the columns of X can combine into, so the residual is orthogonal to each column at once, and X.T @ r holds those dot products, one per column. np.allclose, the array comparison from 2.2.3, checks them all in one call, and its fixed allowance is enough at this size.

import numpy as np

area = np.array([50.0, 60.0, 80.0, 100.0, 120.0])
price = np.array([198.0, 210.0, 285.0, 332.0, 405.0])
X = np.column_stack([np.ones(len(area)), area])
w, ssr, rank, sv = np.linalg.lstsq(X, price)

pred = X @ w                    # 40 + 3 * area, row by row
r = price - pred
print(np.round(pred, 2))        # [190. 220. 280. 340. 400.]
print(np.round(r, 2))           # [  8. -10.   5.  -8.   5.]
print(np.round(r @ r, 6))       # 278.0
print(np.allclose(X.T @ r, 0))  # True

The weights read as a base price of 40 thousand plus 3 thousand per square meter, and the residuals are the gaps no line could close. Why squared residuals are the thing minimized, rather than their absolute sizes, is Module 3's loss lesson, 3.1.3.

DON'T:

• DON'T leave out the column of ones. np.linalg.lstsq(area[:, np.newaxis], price) fits the line through the origin, the projection from the first subtitle, returning the weight [3.44444444] and a sum of squared residuals of [989.11111111] against [278.] with the intercept.
• DON'T unpack lstsq into one name. w = np.linalg.lstsq(X, price) binds the whole tuple, and X @ w then raises ValueError: setting an array element with a sequence. The requested array has an inhomogeneous shape after 1 dimensions. The detected shape was (4,) + inhomogeneous part.
• DON'T read the second value as the residuals. ssr holds one number, their sum of squares, and it comes back as an empty array when the columns are not independent, which 2.3.5 runs, or when the table has no more rows than columns. The residuals are price - X @ w.

The Normal Equations

What. The normal equations, X.T X w = X.T y, are the last subtitle's condition X.T r = 0 written out in the weights: a square system with one equation per weight, which np.linalg.solve answers. Why. They are how least squares is written in textbooks, papers and derivations, and solving them reproduces the weights lstsq found, which shows the right angle is the whole condition. When. When the formula turns up in a derivation or in someone else's code and you need to recognize or reproduce it. To fit a model, lstsq takes X itself.

Write the residual as y - X w inside the orthogonality condition, X.T (y - X w) = 0, multiply out to X.T y - X.T X w = 0, and move one term across: X.T X w = X.T y. X.T @ X is k by k, one row per weight, and X.T @ y has k entries, an ordinary square system. np.linalg.solve(A, c) returns the w for which A @ w equals c, up to rounding.

The textbook writes the answer as w = (X.T X)^-1 X.T y, where ^-1 marks the inverse matrix, np.linalg.inv in NumPy. The inverse of A is the matrix that solves A against every column of the identity, 2.3.2's np.eye, so it solves k systems where one was asked for, and then needs a matrix product on top.

import numpy as np

area = np.array([50.0, 60.0, 80.0, 100.0, 120.0])
price = np.array([198.0, 210.0, 285.0, 332.0, 405.0])
X = np.column_stack([np.ones(len(area)), area])
w, ssr, rank, sv = np.linalg.lstsq(X, price)

XtX = X.T @ X                     # (2, 2), one row per weight
Xty = X.T @ price                 # (2,)
w_ne = np.linalg.solve(XtX, Xty)  # the w where XtX @ w equals Xty, up to rounding
print(np.round(w_ne, 4))          # [40.  3.]
print(np.allclose(w_ne, w))       # True

inv = np.linalg.inv(XtX)
# True, the inverse is XtX solved against each column of the identity
print(np.allclose(inv, np.linalg.solve(XtX, np.eye(2))))
print(np.allclose(inv @ Xty, w))  # True

On five listings and two weights every route agrees. The extra rounding shows on a column with large values beside the column of ones, such as a raw construction year, values near 2000. Putting each answer back into its equations measures how far it misses them, and :.0e (the format spec for scientific notation with no digits after the point) prints each miss.

import numpy as np

price = np.array([198.0, 210.0, 285.0, 332.0, 405.0])
year = np.array([1962.0, 1974.0, 1987.0, 1999.0, 2011.0])
Xy = np.column_stack([np.ones(5), year])
A = Xy.T @ Xy
rhs = Xy.T @ price

w_solve = np.linalg.solve(A, rhs)
w_inv = np.linalg.inv(A) @ rhs
# the largest gap between A @ w and rhs
miss_solve = np.abs(A @ w_solve - rhs).max()
miss_inv = np.abs(A @ w_inv - rhs).max()
print(f"{miss_solve:.0e}")  # 9e-10
print(f"{miss_inv:.0e}")    # 1e-05

The answer from solve lands within a few billionths of the right-hand side, or exactly on it, and the one built on the inverse misses by about a hundred-thousandth, thousands of times farther or more. As with the sliver earlier, your digits for misses this small can differ, and the miss from solve stays the smaller one. How close a matrix like this A is to having no inverse at all is what 2.3.5 measures, with the condition number.

DON'T:

• DON'T type the textbook formula np.linalg.inv(X.T @ X) @ X.T @ y into code. It builds an inverse to use once and misses its own equations by more, as the year column showed, while np.linalg.solve(X.T @ X, X.T @ y) is the same formula done directly.
• DON'T hand solve the design matrix itself. np.linalg.solve(X, price) raises numpy.linalg.LinAlgError: Last 2 dimensions of the array must be square, because five equations in two unknowns is the least squares problem, and solve only takes a square system.
• DON'T drop X.T from the right-hand side. np.linalg.solve(X.T @ X, price) raises ValueError: solve1: Input operand 1 has a mismatch in its core dimension 0, with gufunc signature (m,m),(m)->(m) (size 5 is different from 2), five prices against two equations.

Lab Exercises

Lab Exercise 1: Price per Room Through the Origin

The five listings in this lesson have 2, 3, 3, 4 and 5 rooms. Project their price vector onto the rooms vector to get a price per room for a model with no intercept, then print that price per room and the five fitted prices, both to two decimals.

Compute the residual and test its dot product with the rooms vector twice, once with == and once with np.isclose, and add a comment explaining why the two tests disagree.

Finally, show that your price per room is the closest choice: print the length of the residual for it, for 10 thousand more per room and for 10 thousand less.

Lab Exercise 2: Price on Area and Rooms

Build the design matrix for the same five listings from a column of ones, the areas and the rooms, fit it with np.linalg.lstsq and unpack all four values. Print the weights to four decimals, the rank and the sum of squared residuals.

Recompute the sum of squared residuals from the residual vector and confirm it matches the one lstsq returned, then confirm the residual is orthogonal to every column. Solve the normal equations with np.linalg.solve and confirm they give the same weights.

Finally, predict the price of a listing with 90 square meters and 4 rooms, and add a comment saying what the rooms weight means for two listings of the same area.

import numpy as np

rooms = np.array([2.0, 3.0, 3.0, 4.0, 5.0])
price = np.array([198.0, 210.0, 285.0, 332.0, 405.0])

k = (price @ rooms) / (rooms @ rooms)  # 5234 / 63, a price per room
print(np.round(k, 2))                  # 83.08
print(np.round(k * rooms, 2))          # [166.16 249.24 249.24 332.32 415.4 ]

r = price - k * rooms
print(r @ rooms == 0)                  # False
print(np.isclose(r @ rooms, 0))        # True
# r @ rooms is 2.2737367544323206e-13: 5234 / 63 has no exact float,
# so the residual misses the right angle by rounding, and only == minds

for trial in [k - 10, k, k + 10]:
    gap = np.linalg.norm(price - trial * rooms)
    print(np.round(trial, 2), np.round(gap, 2))

73.08 101.2
83.08 62.77
93.08 101.2

import numpy as np

area = np.array([50.0, 60.0, 80.0, 100.0, 120.0])
rooms = np.array([2.0, 3.0, 3.0, 4.0, 5.0])
price = np.array([198.0, 210.0, 285.0, 332.0, 405.0])

# ones, area, rooms: one row per listing, shape (5, 3)
X = np.column_stack([np.ones(5), area, rooms])
w, ssr, rank, sv = np.linalg.lstsq(X, price)
print(np.round(w, 4))                 # [ 44.1695   3.6407 -16.678 ]
print(rank)                           # 3
print(ssr)                            # [177.93220339]

r = price - X @ w
print(np.isclose(r @ r, ssr[0]))      # True
print(np.allclose(X.T @ r, 0))        # True

w_ne = np.linalg.solve(X.T @ X, X.T @ price)
print(np.allclose(w_ne, w))           # True

listing = np.array([1.0, 90.0, 4.0])  # the 1 meets the intercept
print(np.round(listing @ w, 2))       # 305.12
# -16.678: of two listings with the same area, the one with an extra
# room is predicted 16.678 thousand lower, its rooms being smaller

Theory Assessment

Q1: Doubling the Direction

In the first subtitle, b = [4, 2] gave coef 0.5, proj [2. 1.] and a scalar projection of 2.2361. Replace b with 2 * b, which is [8, 4], and say without running anything what each of the three becomes. Then say which property of b the projection depends on, and what that means for a model's weight when a feature's units change.

Hint: Count how many copies of the longer arrow it takes to reach the same point, and ask whether the line through b moved at all.

Q2: A Sliver Away from Zero

The residual of price projected onto area left r @ area a few trillionths away from zero. Explain why that is not 0.0 although the projection is correct, why the right angle means no other price per square meter could fit better, and what r @ r measures for the five listings.

Hint: For the first part, ask whether 31/9 fits in a float. For the second, remember the sweep along the line, and what the distance did on either side of t = 0.5.

Q3: Four Values and a Column of Ones

Name the four values np.linalg.lstsq(X, price) returned for the five listings and say what each holds. Say which one gets mistaken for the residuals and how the residuals are actually obtained. Then explain what the column of ones changes about the model, using the two sums of squared residuals this lesson printed.

Hint: One of the four is a single number inside an array. For the second part, ask what price each model predicts for a listing with an area of zero.

Q4: Residuals That Sum to Zero

Read the following without running it, then state what the two prints show.

import numpy as np

area = np.array([50.0, 60.0, 80.0, 100.0, 120.0])
price = np.array([198.0, 210.0, 285.0, 332.0, 405.0])

X = np.column_stack([np.ones(5), area])
w, ssr, rank, sv = np.linalg.lstsq(X, price)
r_line = price - X @ w

c = (price @ area) / (area @ area)
r_origin = price - c * area

print(np.isclose(r_line.sum(), 0))
print(np.round(r_origin.sum(), 2))

A teammate reads the first print as proof that the model with the intercept has no error, and the second as proof that the fit through the origin is wrong. Explain what each sum measures, why both conclusions fail, and which number does measure how well each model fits.

Hint: The sum of a vector's entries is its dot product with a vector of ones. Ask which of the two design matrices had a column of ones, and what least squares guarantees about every column it was given.

import numpy as np

a = np.array([1.0, 3.0])
b2 = np.array([8.0, 4.0])  # 2 * b

coef = (a @ b2) / (b2 @ b2)
print(coef)                # 0.25
print(coef * b2)           # [2. 1.]
# 2.2361
print(np.round((a @ b2) / np.linalg.norm(b2), 4))

import numpy as np

area = np.array([50.0, 60.0, 80.0, 100.0, 120.0])
price = np.array([198.0, 210.0, 285.0, 332.0, 405.0])
X = np.column_stack([np.ones(5), area])

shuffled = price[[3, 0, 4, 1, 2]]  # the right prices, wrong listings
w, ssr, rank, sv = np.linalg.lstsq(X, shuffled)
# True
print(np.isclose((shuffled - X @ w).sum(), 0))
print(np.round(ssr, 2))            # [29373.49]

Assignment

Assignment: A Line, a Curve and np.polyfit

np.polyfit(x, y, deg) fits a polynomial of degree deg (a weighted sum of the powers of x up to x ** deg) by least squares and returns its coefficients highest power first, so a straight line comes back as [slope, intercept], its slope being the weight on x. Run it beside np.linalg.lstsq on the lesson's five listings, with their prices in plain currency units as the second subtitle priced them, and reconcile the two.

Show that the straight line from np.polyfit matches the lesson's line, 40 thousand plus 3 thousand per square meter, once reversed. Then fit the curve price = w0 + w1 * area + w2 * area ** 2 by giving lstsq a third column, area ** 2, and print its three weights and its sum of squared residuals beside the line's. Prove the curve's residual orthogonal to all three columns, print what the plain np.allclose(X2.T @ r2, 0) says, and account for the difference. Confirm that np.polyfit(area, price, 2) returns the same three weights.

Finish by pricing a 90 square meter listing with both models, and by explaining why the curve's sum of squared residuals could not have come out larger than the line's, whatever the five prices had been.

Hints:

• A design matrix built as ones, area, area ** 2 holds its weights in the opposite order from np.polyfit, and [::-1] reverses an array.
• One of the curve's three weights, set to zero, turns it back into the line. Ask what least squares does with an option that is always open to it.

import numpy as np

area = np.array([50.0, 60.0, 80.0, 100.0, 120.0])
# the lesson's prices in plain currency units, not thousands
price = np.array([198.0, 210.0, 285.0, 332.0, 405.0]) * 1000
X = np.column_stack([np.ones(5), area])
w, ssr, rank, sv = np.linalg.lstsq(X, price)

line = np.polyfit(area, price, 1)     # slope first, then intercept
print(line)                           # [ 3000. 40000.]
print(np.allclose(line[::-1], w))     # True

X2 = np.column_stack([np.ones(5), area, area ** 2])
w2, ssr2, rank2, sv2 = np.linalg.lstsq(X2, price)
# 82268.04 1899.29 6.5028
print(f"{w2[0]:.2f} {w2[1]:.2f} {w2[2]:.4f}")
print(f"{ssr2[0]:.0f} {ssr[0]:.0f}")  # 225977795 278000000

r2 = price - X2 @ w2
print(np.allclose(X2.T @ r2, 0))      # False
# the cosine between r2 and each of the three columns
cos = (X2.T @ r2) / (np.linalg.norm(X2, axis=0) * np.linalg.norm(r2))
print(np.allclose(cos, 0))            # True

curve = np.polyfit(area, price, 2)
print(np.allclose(curve[::-1], w2))   # True

# 310000.0
print(np.round(np.array([1.0, 90.0]) @ w, 2))
# 305876.29, the curve's row is 1, 90 and 90 ** 2
print(np.round(np.array([1.0, 90.0, 8100.0]) @ w2, 2))

Additional Resources

Optional, and nothing later depends on either. Worth opening if you want the matrix that does the projecting, or the methods solvers use instead of the normal equations.

• YouTube - "16. Projection Matrices and Least Squares" by MIT OpenCourseWare, a lecture by Gilbert Strang. The matrix that projects any vector onto a set of columns, and a line fitted to three points seen two ways, as vertical gaps and as one vector projected in three dimensions.
• Wikipedia - "Numerical methods for linear least squares", section "Orthogonal decomposition methods". How solvers reach the least squares weights without forming X.T X at all.

Summary

One geometric fact carries this lesson: the closest point on a line, or on everything a set of columns can build, is the one whose leftover meets it at a right angle. Projection onto one vector is that fact with one column, least squares is the same fact with one column per weight, and the normal equations are the right angle written out as algebra. Three names, one picture.

The residual is where the picture turns practical. It is the part of the target the features could not reach, and the right angle says it shares nothing with any of them. That same right angle makes some facts about residuals mechanical rather than informative, so residuals that sum to zero only report that a column of ones was there. A fit is judged by the sum of squared residuals it leaves, and even that number speaks only for the listings it was fitted to.

Floats keep the picture honest in their own way. A right angle computed in floating point is a right angle up to rounding, so every check carries a tolerance, and a tolerance has a scale, which is why the cosine survives where a fixed allowance fails. The same care separates solve from the inverse, both correct on paper and the inverse measurably farther from its equations once a real column with large values arrives. What remains open is the case where the columns stop being independent, and the right angle stops picking out a single set of weights.

Cheat Sheet

Part A - Reference Table

Part B - Compact Code Block

import numpy as np

proj = (a @ b) / (b @ b) * b              # closest point on b's line
r = a - proj                              # residual, orthogonal to b
np.isclose(r @ b, 0)                      # never r @ b == 0

X = np.column_stack([np.ones(len(x)), x])
w, ssr, rank, sv = np.linalg.lstsq(X, y)
r = y - X @ w                             # residuals, not ssr
w_ne = np.linalg.solve(X.T @ X, X.T @ y)  # same w if columns independent, no inv