Solving 20,433 equations that have no solution

Predict house value from district features, on real 1990 census data.

d = housing.dropna()                       # 20433 rows remain
X = np.column_stack([np.ones(len(d)), d[feats].to_numpy(float)])   # (20433, 7)
w = np.linalg.pinv(X) @ y
w_lstsq, *_ = np.linalg.lstsq(X, y, rcond=None)
np.allclose(w, w_lstsq)                     # True
rmse = np.sqrt(((X @ w - y) ** 2).mean())   # ≈ 75,980

X is 20,433 × 7 — very tall, so np.linalg.inv cannot even be called. There is no exact solution: no straight line passes through 20,433 points. The pseudoinverse gives the best possible answer instead, and lstsq agrees exactly because it solves the same problem.