Lab U4: Regression as Projection

Unit: Unit 4, Orthogonality and Least Squares
Role: Required
Textbook sections: orthogonal projection; least squares as projection; QR factorization and least squares; applications and computation recap
Core path: projection, design matrices, least-squares coefficients, fitted outputs, residual orthogonality, normal equations, QR, and computational checks
Extension path: adding features and redundant columns as short interpretation checks

Computational tools used in this lab

Before starting, review these parts of Appendix B, NumPy and SymPy Quick Reference for the Labs:

• Appendix B.2: NumPy arrays, vectors, matrices, and shapes
• Appendix B.4: Elementwise arithmetic versus linear algebra
• Appendix B.8: Constructing common arrays and design matrices
• Appendix B.9: Numerical checks and roundoff
• Appendix B.10: NumPy linear algebra commands

The goal is to interpret the mathematical computation, not to memorize every command.

Part 0. Warm-up: shapes, transposes, and @

Math reminder. In NumPy, @ means linear algebra multiplication. With two compatible one-dimensional arrays it is a dot product. With a matrix and a compatible vector it is a matrix-vector product.

Predict before running. Which objects below are vectors? Which object is a matrix? What operation does x @ u perform? Why is A.T @ b defined, and what shape should it have?

import numpy as np

np.set_printoptions(precision=3, suppress=True)

u = np.array([3.0, 4.0])
x = np.array([2.0, 5.0])

A = np.array([
    [1.0, 0.0],
    [0.0, 1.0],
    [1.0, 1.0],
])
b = np.array([1.0, 2.0, 4.0])

u.shape, x.shape, A.shape, b.shape, x @ u, A.T @ b

Run and compare. u, x, and b are vectors; A is a matrix. The expression x @ u is a dot product. Since A is 3 x 2, A.T is 2 x 3, so A.T @ b is defined and produces a 2-vector.

Interpretation check. Reading @ starts with the objects: vector with vector gives a dot product, while matrix with vector gives a matrix-vector product.

Common mistake. A.T @ b is not entry-by-entry multiplication. It computes dot products between the columns of A and the vector b.

Part 1. Projection onto a line

Math reminder. The projection of a vector (\mathbf{x}) onto the line spanned by a nonzero vector (\mathbf{u}) is [ \operatorname{proj}_{\operatorname{span}{\mathbf{u}}}(\mathbf{x})

\frac{\mathbf{x}\cdot\mathbf{u}}{\mathbf{u}\cdot\mathbf{u}}\mathbf{u}. ] The residual is (\mathbf{r}=\mathbf{x}-\operatorname{proj}_{\operatorname{span}{\mathbf{u}}}(\mathbf{x})).

Predict before running. Which variable should store the projection? Which variable should store the residual? Should r @ u be zero, positive, or negative?

u = np.array([3.0, 4.0])
x = np.array([2.0, 5.0])

y = (x @ u) / (u @ u) * u
r = x - y

y, r, r @ u

Run and compare. The variable y stores the projection, and r stores (\mathbf{x}-\mathbf{y}). The value r @ u should be close to zero.

Interpretation check. The residual after projecting onto span{u} is orthogonal to the line direction. That is why r @ u checks the geometry.

Answer sketch. The code represents (((\mathbf{x}\cdot\mathbf{u})/(\mathbf{u}\cdot\mathbf{u}))\mathbf{u}). For these vectors, the projection is approximately [3.12, 4.16] and the residual is approximately [-1.12, 0.84].

Part 2. Orthogonal residuals and A.T @ r

Math reminder. A vector (\mathbf{r}) is orthogonal to every column of (A) exactly when [ A^T\mathbf{r}=\mathbf{0}. ] This is the computational form of (\operatorname{col}(A)^\perp=\operatorname{null}(A^T)).

Predict before running. How many columns does A have? How many entries should A.T @ r have? Does A.T @ r check orthogonality to the rows or to the columns of A?

A = np.array([
    [1.0, 0.0],
    [0.0, 1.0],
    [1.0, 1.0],
])

r = np.array([-1/3, -1/3, 1/3])

A.T @ r

Run and compare. The matrix A has two columns, so A.T @ r has two entries.

Interpretation check. The output checks dot products with the columns of A. If A.T @ r is the zero vector, then (\mathbf{r}\in\operatorname{null}(A^T)=\operatorname{col}(A)^\perp).

Common mistake. The transpose matters: A.T @ r checks column orthogonality, while A @ r is not the right shape here.

Part 3. Unit 2 target revisited

Math reminder. Unit 2 asked whether a target output was reachable. Least squares asks for the closest reachable output when the target is not reachable.

For this matrix, [ A\mathbf{x}=(x_1,x_2,x_1+x_2). ]

Predict before running. Why does (A\mathbf{x}=\mathbf{b}) have no exact solution for b = [1, 2, 4]? Should the residual be the zero vector? Should A.T @ r be close to the zero vector?

A = np.array([
    [1.0, 0.0],
    [0.0, 1.0],
    [1.0, 1.0],
])
b = np.array([1.0, 2.0, 4.0])

xhat = np.linalg.lstsq(A, b, rcond=None)[0]
bhat = A @ xhat
r = b - bhat

xhat, bhat, r, A.T @ r

Run and compare. The target is not reachable because the third coordinate would need to equal the sum of the first two coordinates, but 4 != 1 + 2.

Interpretation check. The vector bhat is the closest reachable output. The residual r = b - bhat is not zero, but A.T @ r is close to zero because the residual is orthogonal to the column space.

Answer sketch. The least-squares coefficient vector is approximately [1.333, 2.333], the fitted output is approximately [1.333, 2.333, 3.667], and the residual is approximately [-0.333, -0.333, 0.333].

Part 4. Design matrices for line fitting

Math reminder. To fit a line (y=c_0+c_1t), the design matrix has one column for the constant feature and one column for the (t)-feature.

Predict before running. What should the first column of X contain? What should the second column contain? Which entries of beta should be (c_0) and (c_1)?

t = np.array([0.0, 1.0, 2.0])
y = np.array([1.0, 2.0, 2.0])

X = np.column_stack([np.ones_like(t), t])

beta = np.linalg.lstsq(X, y, rcond=None)[0]
yhat = X @ beta
r = y - yhat

X, beta, yhat, r

Run and compare. The first column of X is the constant feature. The second column is the t-feature. The vector beta stores the fitted coefficients [c0, c1].

Interpretation check. The vector yhat contains the fitted outputs at the three input values. The vector r contains the residuals (\mathbf{y}-\hat{\mathbf{y}}).

Answer sketch. This is the design-matrix version of fitting (y=c_0+c_1t) to three data points.

Part 5. Normal equations versus lstsq

Math reminder. The normal equations for the design matrix X are [ X^TX\hat{\beta}=X^T\mathbf{y}. ] When the columns of X are independent, X.T @ X is invertible and we can solve this square system.

Predict before running. Why is X.T @ X square? What normal equation is represented by left @ beta = right? Why should beta_normal and beta_lstsq agree?

left = X.T @ X
right = X.T @ y

beta_normal = np.linalg.solve(left, right)
beta_lstsq = np.linalg.lstsq(X, y, rcond=None)[0]

left, right, beta_normal, beta_lstsq

Run and compare. Since X is 3 x 2, X.T @ X is 2 x 2. The equation is (X^TX\hat{\beta}=X^T\mathbf{y}).

Interpretation check. The vectors beta_normal and beta_lstsq agree because they compute the same least-squares coefficients for this independent-column design matrix.

Common mistake. np.linalg.solve(left, right) is appropriate here because left = X.T @ X is square and invertible in this example. That is not guaranteed for every design matrix.

Part 6. Residual orthogonality and roundoff

Math reminder. In exact arithmetic, the least-squares residual satisfies (X^T\mathbf{r}=\mathbf{0}). In floating-point arithmetic, we usually check whether the entries are close to zero.

Predict before running. Does np.allclose(X.T @ r, np.zeros(2)) check that the residual is zero? What does it actually check? If the residual norm is nonzero, did least squares fail?

r = y - X @ beta_lstsq

X.T @ r, np.allclose(X.T @ r, np.zeros(2)), np.linalg.norm(r)

Run and compare. The expression np.allclose(X.T @ r, np.zeros(2)) checks residual orthogonality, not whether the residual itself is zero.

Interpretation check. A nonzero residual norm is normal when the data vector is not exactly in the column space. Least squares has not failed; it found the closest vector in the column space.

Debug. A student writes X @ r to check residual orthogonality. What is wrong?

Answer sketch. The residual lives in output space, and the check should use dot products with the columns of X, encoded as X.T @ r. Also, X @ r has incompatible shapes here.

Part 7. QR least squares

Math reminder. If (X=QR) with Q having orthonormal columns, then least squares can be solved through the triangular system [ R\hat{\beta}=Q^T\mathbf{y}. ]

Predict before running. What should Q.T @ Q be close to? Why does the code solve a system with R instead of a system with X? Should beta_qr and beta_lstsq agree?

Q, R = np.linalg.qr(X)

beta_qr = np.linalg.solve(R, Q.T @ y)
beta_lstsq = np.linalg.lstsq(X, y, rcond=None)[0]

Q.T @ Q, R, beta_qr, beta_lstsq

Run and compare. The matrix Q.T @ Q should be close to the identity matrix because Q has orthonormal columns.

Interpretation check. QR turns the least-squares problem into the triangular solve (R\hat{\beta}=Q^T\mathbf{y}). The vectors beta_qr and beta_lstsq agree because they compute the same least-squares solution.

Common mistake. QR does not prove that (X\beta=\mathbf{y}) has an exact solution. It computes the closest fit when an exact solution is not available.

Part 8. Adding features

Math reminder. Adding a feature adds a column to the design matrix. That gives least squares more possible fitted vectors.

Predict before running. What new column is added to the design matrix? Can the residual norm increase when a column is added? What should X2.T @ r2 check?

X2 = np.column_stack([np.ones_like(t), t, t**2])

beta2 = np.linalg.lstsq(X2, y, rcond=None)[0]
yhat2 = X2 @ beta2
r2 = y - yhat2

np.linalg.norm(r), np.linalg.norm(r2), beta2, X2.T @ r2

Run and compare. The new column is the quadratic feature t**2. Adding a column gives more possible fitted vectors, so the least-squares residual norm cannot increase.

Interpretation check. In this small example there are three points and three independent columns, so the fit can become exact up to roundoff. The expression X2.T @ r2 checks residual orthogonality for the larger column space.

Common mistake. A smaller residual norm does not automatically mean a model is better for every purpose. Here we are only reading the linear algebra of the fit.

Part 9. Redundant columns and nonunique coefficients

Math reminder. Unit 2 null-space directions return in least squares. If X_red @ z = 0, then changing coefficients by a multiple of z does not change the fitted output.

Predict before running. What would X_red @ z = 0 say about the columns of X_red? If the coefficient vector changes by 5*z, should the fitted vector change?

X_red = np.column_stack([np.ones_like(t), t, 1 + t])
z = np.array([1.0, 1.0, -1.0])

beta_red = np.linalg.lstsq(X_red, y, rcond=None)[0]
beta_changed = beta_red + 5*z

X_red @ z, np.linalg.matrix_rank(X_red), X_red @ beta_red, X_red @ beta_changed

Run and compare. The output X_red @ z is the zero vector, so z is a null-space direction for the design matrix. The rank is less than the number of columns, so the columns are dependent.

Interpretation check. The coefficient vector changes, but the fitted vector does not. Coefficients can be nonunique even when the fitted output is fixed.

Answer sketch. Redundant columns create different coefficient descriptions of the same reachable output.

Part 10. Review code-reading bank

Use these as discussion questions or self-review checks. They are code-reading questions, not requests to write new code.

1. Shape check. If X is 3 x 2 and y has shape (3,), what is the shape of X.T @ y?
2. Operation ID. In x @ u, what operation is being computed?
3. Projection formula. What mathematical formula is represented by (x @ u)/(u @ u)*u?
4. Residual. If y = projection and r = x - y, what does r measure?
5. Orthogonality check. Why should r @ u be close to zero after projecting onto span{u}?
6. Column-space check. What does A.T @ r check?
7. Least squares. What does np.linalg.lstsq(A, b, rcond=None)[0] return?
8. Fitted vector. What does A @ xhat represent?
9. Exact solution. If the residual is nonzero, did least squares fail?
10. Normal equations. What equation is represented by A.T @ A @ xhat = A.T @ b?
11. Design matrix. Why does a line fit use a column of ones?
12. QR check. What should Q.T @ Q be close to?
13. QR solve. Why does QR lead to a solve with R?
14. Roundoff. Why do we use np.allclose instead of ==?
15. Debug. What is wrong with checking A @ r instead of A.T @ r?
16. Redundancy. If A @ z = 0, what happens to the fitted vector when (\hat{\mathbf{x}}) is replaced by (\hat{\mathbf{x}}+5\mathbf{z})?
17. Rank. What does rank say about independent columns?
18. Interpretation. Does A.T @ r = 0 mean the residual vector is zero?

Answer sketches

1. X.T @ y has shape (2,).
2. It is a dot product.
3. It is the projection of x onto the line spanned by u.
4. It measures the part of x left over after subtracting the projection.
5. The residual is orthogonal to the line direction.
6. It checks whether r is orthogonal to every column of A.
7. It returns a least-squares coefficient vector.
8. It is the fitted or reachable output produced by xhat.
9. No. A nonzero residual is expected when the target is not exactly reachable.
10. It represents the normal equations (A^TA\hat{\mathbf{x}}=A^T\mathbf{b}).
11. The column of ones represents the constant term.
12. It should be close to an identity matrix.
13. QR rewrites least squares as a triangular solve (R\hat{\mathbf{x}}=Q^T\mathbf{b}).
14. Floating-point arithmetic usually gives tiny roundoff errors.
15. A.T @ r checks dot products with columns of A; A @ r is the wrong operation and may not even have compatible shapes.
16. The coefficient vector changes, but the fitted vector does not.
17. Full column rank means the columns are independent; smaller rank means there is redundancy.
18. No. It means the residual is orthogonal to the column space, not necessarily zero.

Exam-style check

Given r = b - A @ xhat and small values in A.T @ r, identify the geometric condition being checked and explain why this does not prove r = 0.