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, closest reachable outputs, design matrices, least-squares coefficients, fitted vectors, residual orthogonality, normal equations, reduced QR, and numerical checks
Assessment note: This notebook supplements the lecture notes. Predict each output, run the cell, and interpret both the mathematics and the NumPy syntax.

Submission note: No code submission is expected; this lab supports in-class activities and guided review.

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 read short computations as mathematical notation and to understand what each NumPy command returns.

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 expression computes the scalar projection coefficient? What array should xproj store? What should r store? Should r @ u be exactly zero in floating-point arithmetic?

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

coefficient = (x @ u) / (u @ u)
xproj = coefficient * u
r = x - xproj

xproj, r, r @ u, np.allclose(r @ u, 0.0)

Run and compare. The scalar coefficient is

[ \frac{\mathbf{x}\cdot\mathbf{u}} {\mathbf{u}\cdot\mathbf{u}} = \frac{26}{25}. ]

The projection and residual are

[ \widehat{\mathbf{x}} = \begin{bmatrix} 3.12\\ 4.16 \end{bmatrix}, \qquad \mathbf{r} = \begin{bmatrix} -1.12\\ 0.84 \end{bmatrix}. ]

Syntax check. The expression x @ u is a dot product. Multiplying the scalar coefficient by u produces the projected vector.

Interpretation check. In exact arithmetic, (\mathbf{r}\cdot\mathbf{u}=0). Floating-point arithmetic may return a tiny number such as -8.88e-16, while np.allclose returns True.

Common mistake. The residual is x - xproj, not xproj - x.

Part 2. A.T checks column orthogonality

Math reminder. “A^T checks column orthogonality” says that

[ A^T\mathbf{r} = \begin{bmatrix} \mathbf{a}_1^T\mathbf{r}\\ \vdots\\ \mathbf{a}_n^T\mathbf{r} \end{bmatrix}. ]

Thus

[ A^T\mathbf{r}=\mathbf{0} ]

exactly when (\mathbf{r}) is orthogonal to every column of (A).

Predict before running. How many columns does A have? What should A.shape[1] return? What shape should A.T @ r have? Why does the zero array in the numerical check need the same length?

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

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

orthogonality = A.T @ r

orthogonality, np.allclose(
    orthogonality,
    np.zeros(A.shape[1]),
)

Run and compare. The matrix A has two columns, so A.shape[1] is 2 and A.T @ r has two entries.

Syntax check.

np.zeros(A.shape[1])

constructs a zero vector with one entry for each column of A.

Interpretation check. The output checks the two dot products between the residual and the columns of (A). The Boolean result is True, so

[ \mathbf{r} \in \operatorname{null}(A^T) = \operatorname{col}(A)^\perp. ]

Common mistake. The transpose matters. A.T @ r collects dot products with the columns of A; A @ r is not the correct operation.

Part 3. Closest reachable output

Math reminder. “Unit 2 target revisited” begins with a target (\mathbf{b}) that is not in (\operatorname{col}(A)). Least squares finds a coefficient vector (\widehat{\mathbf{x}}) whose fitted vector

[ \widehat{\mathbf{b}} = A\widehat{\mathbf{x}} ]

is the closest reachable output. The residual is

[ \mathbf{r} = \mathbf{b}-\widehat{\mathbf{b}}. ]

For this matrix,

[ A\mathbf{x} = \begin{bmatrix} x_1\\ x_2\\ x_1+x_2 \end{bmatrix}. ]

Predict before running. Why is b = [1, 2, 4] not reachable? What does [0] select from the value returned by np.linalg.lstsq? Which variable should store the coefficient vector? Which should store the closest reachable output? Should the residual be zero?

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, np.linalg.norm(r)**2

Run and compare. The first two target coordinates would force (x_1=1) and (x_2=2), but then the third output coordinate would be (3), not (4).

The least-squares coefficient vector is

[ \widehat{\mathbf{x}} = \begin{bmatrix} 4/3\\ 7/3 \end{bmatrix}. ]

The fitted vector and residual are

[ \widehat{\mathbf{b}} = \begin{bmatrix} 4/3\\ 7/3\\ 11/3 \end{bmatrix}, \qquad \mathbf{r} = \begin{bmatrix} -1/3\\ -1/3\\ 1/3 \end{bmatrix}. ]

Also,

[ A^T\mathbf{r}=\mathbf{0}, \qquad |\mathbf{r}|^2=\frac13. ]

Syntax check. The function np.linalg.lstsq returns several objects. The index [0] selects the least-squares coefficient vector.

Interpretation check. xhat is a coefficient vector in (\mathbb{R}^2). bhat is the closest reachable output in (\mathbb{R}^3). The residual is nonzero because the original target is not reachable.

Part 4. A line model as a design matrix

Math reminder. “A line model as a matrix product” writes the predictions

[ \widehat y=c_0+c_1t ]

as a matrix-vector product. The design matrix contains one column for the constant feature and one column for the (t)-feature.

The figure “Two views of the same line fit” distinguishes vertical residuals in the (t)-(y) data plot from orthogonal projection in output space.

The activity “Fit a line to three data points” carries out this line fit by hand.

Predict before running. What array should np.ones_like(t) produce? What does np.column_stack do? What should the shape of X be? Why should c have two entries? Which variables should store the fitted vector and residual?

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])

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

X, c, yhat, r, X.T @ r, np.linalg.norm(r)**2

Run and compare. The design matrix is

[ X = \begin{bmatrix} 1&0\\ 1&1\\ 1&2 \end{bmatrix}. ]

The coefficient vector is

[ \widehat{\mathbf{c}} = \begin{bmatrix} 7/6\\ 1/2 \end{bmatrix}, ]

so the fitted line is

[ \widehat y = \frac76+\frac12t. ]

The fitted vector and residual are

[ \widehat{\mathbf{y}} = \begin{bmatrix} 7/6\\ 5/3\\ 13/6 \end{bmatrix}, \qquad \mathbf{r} = \begin{bmatrix} -1/6\\ 1/3\\ -1/6 \end{bmatrix}. ]

Also,

[ X^T\mathbf{r}=\mathbf{0}, \qquad |\mathbf{r}|^2=\frac16. ]

Syntax check. np.ones_like(t) produces a length-three array of ones. np.column_stack uses that array and t as the two columns of X.

Interpretation check. c stores the intercept and slope. yhat stores one fitted output for each input value. In the data plot, the residuals are vertical differences (y_i-\widehat y_i).

Part 5. Normal equations and lstsq

Math reminder. “Least squares and the normal equations” says that the least-squares coefficients satisfy

[ X^TX\widehat{\mathbf{c}} = X^T\mathbf{y}. ]

When the columns of (X) are linearly independent, (X^TX) is invertible and this square system has a unique solution.

Predict before running. Why is X.T @ X square? What mathematical objects are stored in normal_matrix and normal_rhs? Why is np.linalg.solve appropriate here? Why should c_normal and c_lstsq agree?

normal_matrix = X.T @ X
normal_rhs = X.T @ y

c_normal = np.linalg.solve(normal_matrix, normal_rhs)
c_lstsq = np.linalg.lstsq(X, y, rcond=None)[0]

(
    normal_matrix,
    normal_rhs,
    c_normal,
    c_lstsq,
    np.allclose(c_normal, c_lstsq),
)

Run and compare.

[ X^TX = \begin{bmatrix} 3&3\\ 3&5 \end{bmatrix}, \qquad X^T\mathbf{y} = \begin{bmatrix} 5\\ 6 \end{bmatrix}. ]

Both methods return

[ \widehat{\mathbf{c}} = \begin{bmatrix} 7/6\\ 1/2 \end{bmatrix}, ]

and the final Boolean comparison is True.

Syntax check. np.linalg.solve(normal_matrix, normal_rhs) solves a square linear system. It is preferable to explicitly constructing an inverse matrix.

Interpretation check. The normal equations and np.linalg.lstsq describe the same least-squares problem. The square solve is valid here because the two columns of (X) are linearly independent.

Common mistake. A general design matrix need not have independent columns. In that case, (X^TX) is not invertible even though least-squares fitted vectors still exist.

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, the entries may be tiny rather than exactly zero.

Predict before running. What does X.shape[1] return? What array does np.zeros(X.shape[1]) construct? What does np.allclose check? Does a nonzero residual norm mean that least squares failed?

r = y - X @ c_lstsq
orthogonality = X.T @ r

(
    orthogonality,
    np.allclose(
        orthogonality,
        np.zeros(X.shape[1]),
    ),
    np.linalg.norm(r),
    np.linalg.norm(r)**2,
)

Run and compare. The orthogonality vector is numerically close to

array([0., 0.])

and the Boolean check is True.

The residual norm and squared residual norm are

[ |\mathbf{r}| = \frac{1}{\sqrt6} \approx 0.408248, \qquad |\mathbf{r}|^2 = \frac16. ]

Syntax check. X.shape[1] is the number of columns of (X), so

np.zeros(X.shape[1])

constructs a zero vector with one entry for each residual-orthogonality condition.

Interpretation check. np.allclose checks whether (X^T\mathbf{r}) is numerically close to zero. It does not check whether (\mathbf{r}) itself is zero.

Debug. X @ r is not the residual-orthogonality check. It has incompatible shapes here, and it does not compute dot products with the columns of (X).

Part 7. Least squares with reduced QR

Math reminder. QR reuses the orthonormal directions constructed by the “Gram–Schmidt algorithm.” The algorithm “Least squares with QR” solves

[ R\widehat{\mathbf{c}} = Q^T\mathbf{y} ]

after computing a reduced factorization

[ X=QR. ]

Practice the projection, subtraction, measurement, and normalization sequence in “Perfectly Normal.”

Predict before running. Why does the code request mode="reduced"? What should the shapes of Q and R be? What do the two structural Boolean checks test? Why does the code solve a system with R? Should c_qr and c_lstsq agree?

Q, R = np.linalg.qr(X, mode="reduced")

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

(
    Q.shape,
    R.shape,
    np.allclose(Q.T @ Q, np.eye(Q.shape[1])),
    np.allclose(Q @ R, X),
    c_qr,
    c_lstsq,
    np.allclose(c_qr, c_lstsq),
)

Run and compare. Since (X) has shape (3\times2), reduced QR gives

Q.shape == (3, 2)
R.shape == (2, 2)

The first structural check verifies

[ Q^TQ=I_2, ]

and the second verifies

[ QR=X. ]

Both coefficient computations return values close to

array([1.16666667, 0.5])

and the final comparison is True.

Syntax check. The assignment

Q, R = np.linalg.qr(X, mode="reduced")

stores the two arrays returned by the QR routine. The expression

np.eye(Q.shape[1])

constructs the (2\times2) identity matrix.

Interpretation check. QR replaces the rectangular least-squares problem by the square upper-triangular system

[ R\widehat{\mathbf{c}} = Q^T\mathbf{y}. ]

Agreement with np.linalg.lstsq does not mean that (X\mathbf{c}=\mathbf{y}) has an exact solution.

Common mistake. Numerical software may return signs in (Q) and (R) that differ from a hand Gram–Schmidt calculation. The checks

[ Q^TQ=I \qquad\text{and}\qquad QR=X ]

are the relevant tests.

Part 8. Adding a feature enlarges the fitted space

Math reminder. “A line model as a matrix product” uses the two feature columns (1) and (t). Adding the column (t^2) changes the model from

[ \widehat y=c_0+c_1t ]

to

[ \widehat y=c_0+c_1t+c_2t^2. ]

The original line fits are still available by choosing (c_2=0). Adding a feature therefore cannot increase the smallest possible squared residual.

Predict before running.

1. What array does t**2 produce?

2. What does

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

   construct?

3. What should the shapes of X_line and X_quad be?

4. What should their ranks be for these three input values?

5. Can the quadratic model fit all three observed outputs exactly?

6. What should X_quad.T @ r_quad check?

7. Why can the minimum squared residual not increase when the new feature is added?

8. Does an exact fit to these three observations prove that the quadratic model will predict new outputs better?

X_line = X

c_line = np.linalg.lstsq(X_line, y, rcond=None)[0]
yhat_line = X_line @ c_line
r_line = y - yhat_line

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

c_quad = np.linalg.lstsq(X_quad, y, rcond=None)[0]
yhat_quad = X_quad @ c_quad
r_quad = y - yhat_quad

(
    X_line.shape,
    X_quad.shape,
    np.linalg.matrix_rank(X_line),
    np.linalg.matrix_rank(X_quad),
    c_quad,
    yhat_quad,
    np.linalg.norm(r_line)**2,
    np.linalg.norm(r_quad)**2,
    np.allclose(
        X_quad.T @ r_quad,
        np.zeros(X_quad.shape[1]),
    ),
)

Run and compare. The two design matrices have shapes

X_line.shape == (3, 2)
X_quad.shape == (3, 3)

and ranks

np.linalg.matrix_rank(X_line) == 2
np.linalg.matrix_rank(X_quad) == 3

The quadratic design matrix is

[ X_{\mathrm{quad}} = \begin{bmatrix} 1&0&0\\ 1&1&1\\ 1&2&4 \end{bmatrix}. ]

The fitted coefficient vector is

[ \widehat{\mathbf{c}}_{\mathrm{quad}} = \begin{bmatrix} 1\\ 3/2\\ -1/2 \end{bmatrix}, ]

so the fitted quadratic is

[ \widehat y = 1+\frac32t-\frac12t^2. ]

At the three observed inputs,

[ \widehat{\mathbf{y}}_{\mathrm{quad}} = \begin{bmatrix} 1\\ 2\\ 2 \end{bmatrix}. ]

Thus the quadratic residual is numerically close to the zero vector. The squared residual decreases from

[ |\mathbf{r}_{\mathrm{line}}|^2 = \frac16 ]

to a value close to (0).

Syntax check. The expression t**2 squares each entry of t. np.column_stack places the three feature arrays into the columns of one matrix. np.linalg.matrix_rank reports the number of independent columns.

Interpretation check. Adding a feature cannot increase the smallest squared residual because every old line fit remains available by choosing the new quadratic coefficient to be (0). In this example, three independent columns span all of (\mathbb{R}^3), so the three observed outputs can be fit exactly.

The check

X_quad.T @ r_quad

tests residual orthogonality for the larger column space.

Common mistake. An exact fit to these three observations does not by itself show that the quadratic model will predict new outputs better. This computation compares only the residuals at the displayed data points.

Part 9. Redundant features and nonunique coefficients

Math reminder. A design matrix can contain different columns that repeat the same information. If

[ X\mathbf{z}=\mathbf{0}, ]

then

[ X(\widehat{\mathbf{c}}+\mathbf{z}) = X\widehat{\mathbf{c}}. ]

The coefficient vector changes, but the fitted vector does not.

The matrix below uses three feature columns:

[ 1, \qquad t, \qquad 1+t. ]

The third feature is the sum of the first two.

Predict before running.

1. What matrix is constructed as X_red?

2. Why are its columns dependent?

3. What should X_red @ z return for

   z = np.array([1.0, 1.0, -1.0])
   

   ?

4. What should the rank of X_red be?

5. Should c_red and c_changed be the same vector?

6. Should yhat_red and yhat_changed be the same vector?

7. Which object can be nonunique: the fitted vector or the coefficient vector?

8. Should the residual-orthogonality check still hold when the columns are dependent?

X_red = np.column_stack([
    np.ones_like(t),
    t,
    1 + t,
])

z = np.array([1.0, 1.0, -1.0])

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

yhat_red = X_red @ c_red
yhat_changed = X_red @ c_changed
r_red = y - yhat_red

(
    X_red,
    X_red @ z,
    np.linalg.matrix_rank(X_red),
    c_red,
    c_changed,
    yhat_red,
    yhat_changed,
    np.allclose(yhat_red, yhat_changed),
    X_red.T @ r_red,
    np.linalg.norm(r_red)**2,
)

Run and compare. The redundant design matrix is

[ X_{\mathrm{red}} = \begin{bmatrix} 1&0&1\\ 1&1&2\\ 1&2&3 \end{bmatrix}. ]

Its third column is the sum of its first two columns. Therefore

[ X_{\mathrm{red}} \begin{bmatrix} 1\\ 1\\ -1 \end{bmatrix} = \mathbf{0}, ]

and

np.linalg.matrix_rank(X_red) == 2

even though X_red has three columns.

The vectors c_red and c_changed are different, but

np.allclose(yhat_red, yhat_changed)

returns True. Both produce the fitted vector

[ \widehat{\mathbf{y}} = \begin{bmatrix} 7/6\\ 5/3\\ 13/6 \end{bmatrix}. ]

Their residual is the same line-fitting residual, with

[ |\mathbf{r}_{\mathrm{red}}|^2 = \frac16. ]

Syntax check. The expression

c_changed = c_red + 5*z

moves the coefficient vector in a null-space direction. Matrix multiplication then checks whether that coefficient change affects the fitted output.

Interpretation check. The fitted vector is unique, but the coefficient description is not. Any coefficient vector of the form

[ \widehat{\mathbf{c}}+s\mathbf{z} ]

produces the same fitted vector because (X_{\mathrm{red}}\mathbf{z}=\mathbf{0}).

The residual still satisfies

[ X_{\mathrm{red}}^T\mathbf{r}_{\mathrm{red}} = \mathbf{0}. ]

Common mistake. More columns do not necessarily mean more independent features. A column that is a linear combination of earlier columns does not enlarge the column space.

Part 10. Review activities

Activity A. Read a least-squares pipeline

Consider the code pattern

c = np.linalg.lstsq(X, y, rcond=None)[0]
yhat = X @ c
r = y - yhat
orthogonality = X.T @ r
squared_error = np.linalg.norm(r)**2

1. Which variable stores the coefficient vector?

2. Which variable stores the fitted vector?

3. Which variable stores the residual?

4. If (X) is (m\times n), what are the shapes of c, yhat, r, and orthogonality?

5. What does [0] select from the value returned by np.linalg.lstsq?

6. What geometric condition is checked by X.T @ r?

7. Why can r be nonzero even when X.T @ r is numerically zero?

8. What mathematical quantity is computed by np.linalg.norm(r)**2?

Activity B. Diagnose the claims

For each claim, decide whether it is correct. If it is not correct, replace it with a correct statement.

1. Adding a feature column can make the minimum least-squares squared residual larger.

2. If X @ z is the zero vector, then X @ c and X @ (c + z) are the same fitted vector.

3. If Q.T @ Q is the identity matrix, then Q @ Q.T must also be the identity matrix, even when Q is rectangular.

4. A numerical QR factorization must use exactly the same signs as a hand Gram–Schmidt calculation.

5. The code

   np.linalg.solve(R, Q.T @ y)
   

   solves the reduced system

   [ R\widehat{\mathbf{c}} = Q^T\mathbf{y}. ]

6. In one sentence, connect regression, projection, and QR.

Answer sketches

Activity A.

1. c stores the coefficient vector.

2. yhat stores the fitted vector.

3. r stores the residual.

4. If (X) is (m\times n), then c has shape (n,), while yhat and r have shape (m,). The vector orthogonality has shape (n,).

5. The index [0] selects the least-squares coefficient vector.

6. X.T @ r checks whether the residual is orthogonal to every column of (X).

7. A nonzero residual is expected when the target vector is not in (\operatorname{col}(X)). Least squares requires the residual to be orthogonal to the column space, not equal to zero.

8. np.linalg.norm(r)**2 computes the sum of squared residuals.

Activity B.

1. The claim is not correct. Adding a feature cannot increase the minimum least-squares squared residual because the old fitted vectors remain available.

2. The claim is correct. A null-space change in the coefficients does not change the fitted vector.

3. The claim is not correct. A rectangular matrix with orthonormal columns satisfies

   [ Q^TQ=I, ]

   but generally

   [ QQ^T\neq I. ]

4. The claim is not correct. Simultaneous sign changes in one column of (Q) and the corresponding row of (R) leave (QR) unchanged.

5. The claim is correct.

6. Linear regression projects the data vector onto the column space of a design matrix, and QR supplies orthonormal coordinates that reduce the least-squares computation to an upper-triangular solve.