Lab U6: Constraints, PCA, and SVD

Unit: Unit 6, Constrained optimization
Role: Required
Textbook sections: Lagrange multipliers; quadratic forms and maximum stretch; principal component analysis; singular value decomposition; PCA from the SVD; applications and computation recap
Core path: check one- and several-constraint Lagrange conditions; connect quadratic forms with maximum stretch; center data and form covariance; compute PCA scores, reconstructions, residuals, and principal coordinates; compare PCA with regression; connect PCA with SVD; and read rank and the four fundamental subspaces from an SVD
Submission note: No code submission is expected; this lab supports in-class interpretation questions and guided review.

This lab is about reading short NumPy computations as mathematics. The goal is not to write long programs. Predict the output or shape first, run the cell, and then explain what the computation says.

Parts 0–7 are the core path. Part 8 contains review and sanity 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.9: Numerical checks and roundoff
• Appendix B.10: NumPy linear algebra commands
• Appendix B.14: Optional plotting

The goal is to interpret the mathematical computation, not to memorize every command. The plotting command used below only draws the two fitted lines; it does not perform the PCA or regression calculation.

Part 0. Lagrange conditions in code

Math reminder. With one equality constraint, a constrained critical point satisfies [ \nabla f(\mathbf a)=\lambda\nabla g(\mathbf a), \qquad g(\mathbf a)=k. ]

With two equality constraints, the corresponding condition is [ \nabla f(\mathbf a) = \lambda\nabla g_1(\mathbf a)

\mu\nabla g_2(\mathbf a). ]

The code below verifies these equations at candidates already found by a Lagrange-multiplier calculation. It does not replace the step of finding the candidates.

Predict before running.

1. For (f(x,y)=3x+4y) on (x^2+y^2=25), what signs should the multipliers have at ((3,4)) and ((-3,-4))?
2. For (f(x,y,z)=x-y) subject to [ x^2+y^2+z^2=1, \qquad x+y+z=0, ] must (\nabla f) be parallel to each constraint gradient?

import numpy as np

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

# One constraint: f(x, y) = 3x + 4y on x^2 + y^2 = 25.
a_max = np.array([3.0, 4.0])
a_min = -a_max

grad_f_one = np.array([3.0, 4.0])
grad_g_max = 2 * a_max
grad_g_min = 2 * a_min

lambda_max = (grad_f_one @ grad_g_max) / (grad_g_max @ grad_g_max)
lambda_min = (grad_f_one @ grad_g_min) / (grad_g_min @ grad_g_min)

one_max_ok = np.allclose(grad_f_one, lambda_max * grad_g_max)
one_min_ok = np.allclose(grad_f_one, lambda_min * grad_g_min)

# Two constraints: f(x, y, z) = x - y on the unit sphere
# and the plane x + y + z = 0.
a_two = np.array([1.0, -1.0, 0.0]) / np.sqrt(2)

grad_f_two = np.array([1.0, -1.0, 0.0])
grad_g1 = 2 * a_two
grad_g2 = np.ones(3)

# The columns of G are the two constraint gradients.
G = np.column_stack([grad_g1, grad_g2])
multipliers = np.linalg.lstsq(G, grad_f_two, rcond=None)[0]
several_ok = np.allclose(grad_f_two, G @ multipliers)

(
    (a_max @ a_max, grad_f_one @ a_max, lambda_max, one_max_ok),
    (a_min @ a_min, grad_f_one @ a_min, lambda_min, one_min_ok),
    (
        a_two @ a_two,
        a_two.sum(),
        np.linalg.matrix_rank(G),
        multipliers,
        several_ok,
    ),
)

Run and compare. For the one-constraint problem, both candidates satisfy [ x^2+y^2=25. ] Their objective values are (25) and (-25). The multipliers are [ \lambda_{\max}=\frac12, \qquad \lambda_{\min}=-\frac12, ] and both parallel-gradient checks return True.

For the two-constraint problem, the candidate has unit length and coordinate sum (0). The two constraint gradients are independent because np.linalg.matrix_rank(G) returns 2. The computed multipliers are approximately [ \begin{bmatrix} 1/\sqrt2\\ 0 \end{bmatrix}, ] and the linear-combination check returns True.

Interpretation check. With one constraint, the objective gradient is parallel to one constraint gradient. With several constraints, the objective gradient lies in the span of the constraint gradients. It need not be parallel to each one.

Common mistake. The command np.linalg.lstsq is used here only to find the coefficients in a linear combination of two known constraint gradients. It is not solving the constrained optimization problem.

Part 1. Quadratic forms and maximum stretch

Math reminder. For a matrix (A), [ |A\mathbf x|^2 = \mathbf x^TA^TA\mathbf x. ] Thus maximum stretch on the unit circle is a constrained quadratic-form problem for (A^TA).

Consider [ A= \begin{bmatrix} 1&2\\ 2&1 \end{bmatrix}. ]

Predict before running.

1. Is (A^TA) diagonal?
2. Should the most-stretched direction be a coordinate direction?
3. What should the largest singular value be if the largest eigenvalue of (A^TA) is (9)?

A_stretch = np.array([[1.0, 2.0],
                      [2.0, 1.0]])

B = A_stretch.T @ A_stretch

eigvals_stretch, Q_stretch = np.linalg.eigh(B)
order = np.argsort(eigvals_stretch)[::-1]
eigvals_stretch = eigvals_stretch[order]
Q_stretch = Q_stretch[:, order]

singular_values_stretch = np.sqrt(eigvals_stretch)

# Sample many unit inputs as a numerical check.
theta = np.linspace(0.0, 2 * np.pi, 720, endpoint=False)
unit_inputs = np.column_stack([np.cos(theta), np.sin(theta)])
sampled_stretches = np.linalg.norm(unit_inputs @ A_stretch.T, axis=1)

j_max = np.argmax(sampled_stretches)
sampled_direction = unit_inputs[j_max]
sample_alignment = abs(sampled_direction @ Q_stretch[:, 0])

(
    B,
    eigvals_stretch,
    singular_values_stretch,
    sampled_stretches[j_max],
    sampled_direction,
    sample_alignment,
)

Run and compare. The matrix is [ A^TA = \begin{bmatrix} 5&4\\ 4&5 \end{bmatrix}. ] Its eigenvalues, listed from largest to smallest, are (9) and (1). The singular values are therefore (3) and (1).

The sampled maximum stretch is (3), up to roundoff. The absolute dot product between the sampled maximizing direction and the largest-eigenvalue direction is (1), up to roundoff. The two directions therefore agree up to sign.

Interpretation check. The constrained quadratic form [ \mathbf x^TA^TA\mathbf x ] has maximum value (9) on the unit circle. The corresponding matrix stretch is its square root, (3).

Common mistake. The sampling calculation supports the geometric picture but does not prove the exact maximum. The eigenvalue calculation gives the exact result.

Part 2. Centering and covariance

The rows of [ X= \begin{bmatrix} 4&3\\ 3&4\\ 0&1\\ 1&0 \end{bmatrix} ] are four data points.

Math reminder. PCA begins by centering the data: [ \mathbf z_i=\mathbf x_i-\bar{\mathbf x}. ] If the centered vectors are the rows of (Z), then [ C=\frac1nZ^TZ ] is the covariance matrix used in this unit.

Predict before running.

1. What is the shape of xbar?
2. What should Z.mean(axis=0) return?
3. What is the shape of (C)?

X = np.array([[4.0, 3.0],
              [3.0, 4.0],
              [0.0, 1.0],
              [1.0, 0.0]])

xbar = X.mean(axis=0)
Z = X - xbar
C = Z.T @ Z / len(X)

X.shape, xbar, Z, Z.mean(axis=0), C

Run and compare. The mean is [ \bar{\mathbf x} = \begin{bmatrix} 2\\ 2 \end{bmatrix}, ] and [ Z= \begin{bmatrix} 2&1\\ 1&2\\ -2&-1\\ -1&-2 \end{bmatrix}. ] The two feature means of (Z) are both (0), up to roundoff. The covariance matrix is [ C= \begin{bmatrix} 5/2&2\\ 2&5/2 \end{bmatrix}. ]

Interpretation check. The diagonal entries record average squared centered coordinates. The off-diagonal entry records the average product of the two centered coordinates.

Common mistake. Forming X.T @ X / n without first centering generally analyzes direction from the origin, not variation around the data mean.

Part 3. Scores, reconstructions, and residuals

Math reminder. For a unit direction (\mathbf v), [ \mathbf t=Z\mathbf v ] contains the scores, [ \widehat Z ] contains the projected centered data, and [ R=Z-\widehat Z ] contains the reconstruction residuals.

For each row, [ |\mathbf z_i|^2 = t_i^2+|\mathbf r_i|^2. ]

Compare the two orthogonal directions [ \mathbf v_{\mathrm{diag}} = \frac1{\sqrt2} \begin{bmatrix} 1\\ 1 \end{bmatrix}, \qquad \mathbf v_{\mathrm{cross}} = \frac1{\sqrt2} \begin{bmatrix} 1\\ -1 \end{bmatrix}. ]

Predict before running. Which direction should capture more of this data set?

def project_rows_onto_direction(Z, v):
    v = v / np.linalg.norm(v)

    scores = Z @ v
    Zhat = np.outer(scores, v)
    R = Z - Zhat

    captured = np.sum(scores**2)
    missed = np.sum(R**2)
    total = np.sum(Z**2)

    return scores, Zhat, R, captured, missed, total


v_diag = np.array([1.0, 1.0]) / np.sqrt(2)
v_cross = np.array([1.0, -1.0]) / np.sqrt(2)

scores_diag, Zhat_diag, R_diag, captured_diag, missed_diag, total = (
    project_rows_onto_direction(Z, v_diag)
)

scores_cross, Zhat_cross, R_cross, captured_cross, missed_cross, _ = (
    project_rows_onto_direction(Z, v_cross)
)

residuals_ok = np.allclose(R_diag @ v_diag, np.zeros(len(X)))
pythagoras_ok = np.allclose(total, captured_diag + missed_diag)

(
    (captured_diag, missed_diag),
    (captured_cross, missed_cross),
    total,
    scores_diag,
    Zhat_diag,
    R_diag @ v_diag,
    residuals_ok,
    pythagoras_ok,
)

Run and compare. For (\mathbf v_{\mathrm{diag}}), the total squared score is (18) and the total squared reconstruction residual is (2). For (\mathbf v_{\mathrm{cross}}), these values are (2) and (18). The total centered squared length is (20) in either case.

The matrix [ \widehat Z= \begin{bmatrix} 3/2&3/2\\ 3/2&3/2\\ -3/2&-3/2\\ -3/2&-3/2 \end{bmatrix} ] contains the projections onto the diagonal direction. The output R_diag @ v_diag is a list of one residual dot product for each data point; all are numerically zero.

Interpretation check. The fixed total [ 20=18+2 ] shows why maximizing squared scores is equivalent to minimizing squared reconstruction residuals.

Common mistake. scores_diag is a vector of scalar coordinates. Zhat_diag is a matrix of reconstructed centered data vectors. They are not the same mathematical object.

Part 4. First and second principal directions

Math reminder. The first principal direction is a unit eigenvector of (C) for its largest eigenvalue. The second principal direction is an orthogonal eigenvector for the next eigenvalue.

If [ V= \begin{bmatrix} \mathbf v_1&\mathbf v_2 \end{bmatrix}, ] then [ T=ZV ] contains the principal-component scores and [ V^TCV ] is diagonal.

Predict before running.

1. In what order does np.linalg.eigh return the eigenvalues?
2. Which axis should be the first principal direction?
3. What fraction of the total variation should it capture?

eigvals_pca, V_pca = np.linalg.eigh(C)

order = np.argsort(eigvals_pca)[::-1]
eigvals_pca = eigvals_pca[order]
V_pca = V_pca[:, order]

T = Z @ V_pca
cov_in_principal_coordinates = V_pca.T @ C @ V_pca
score_covariance = T.T @ T / len(X)

fraction_first = eigvals_pca[0] / eigvals_pca.sum()

orthogonal_ok = np.allclose(V_pca.T @ V_pca, np.eye(2))
diagonal_ok = np.allclose(
    cov_in_principal_coordinates,
    np.diag(eigvals_pca),
)
score_cov_ok = np.allclose(
    score_covariance,
    cov_in_principal_coordinates,
)

(
    eigvals_pca,
    V_pca,
    T,
    cov_in_principal_coordinates,
    score_covariance,
    fraction_first,
    orthogonal_ok,
    diagonal_ok,
    score_cov_ok,
)

Run and compare. After reversing the order returned by np.linalg.eigh, the covariance eigenvalues are [ \lambda_1=\frac92, \qquad \lambda_2=\frac12. ] The columns of (V) are the first and second principal directions, up to independent sign changes.

Both [ V^TCV ] and [ \frac1nT^TT ] are numerically [ \begin{bmatrix} 9/2&0\\ 0&1/2 \end{bmatrix}. ] The first principal direction captures [ \frac{9/2}{9/2+1/2} = \frac9{10} ] of the total centered variation.

Interpretation check. The two columns of (T) are score vectors. Their average squared entries are (9/2) and (1/2).

Common mistake. Replacing a principal direction by its negative reverses the corresponding scores but does not change the principal axis or the projected data.

Part 5. PCA and linear regression

PCA and linear regression can fit different lines to the same data.

• PCA treats the two feature coordinates symmetrically and minimizes perpendicular reconstruction residuals.
• Regression of (y) on (x) minimizes vertical prediction residuals.

Both fitted lines pass through the sample mean in this example.

Predict before running.

1. What is the slope of the PCA line?
2. Which line should have the smaller vertical sum of squares?
3. Which line should have the smaller perpendicular sum of squares?

import matplotlib.pyplot as plt

v1 = V_pca[:, 0]
pca_slope = v1[1] / v1[0]

# Regression of y on x, including an intercept.
design = np.column_stack([np.ones(len(X)), X[:, 0]])
beta = np.linalg.lstsq(design, X[:, 1], rcond=None)[0]
yhat_regression = design @ beta

# Vertical predictions from the PCA line through the sample mean.
yhat_pca = xbar[1] + pca_slope * (X[:, 0] - xbar[0])

vertical_sse_regression = np.sum((X[:, 1] - yhat_regression)**2)
vertical_sse_pca = np.sum((X[:, 1] - yhat_pca)**2)

# Perpendicular residuals to the regression line.
regression_direction = np.array([1.0, beta[1]])
regression_direction /= np.linalg.norm(regression_direction)

_, _, R_regression_axis, _, perpendicular_sse_regression, _ = (
    project_rows_onto_direction(Z, regression_direction)
)

perpendicular_sse_pca = np.sum(R_diag**2)

print("regression coefficients [b, m]:", beta)
print("PCA slope:", pca_slope)
print("vertical SSE: regression, PCA:",
      vertical_sse_regression, vertical_sse_pca)
print("perpendicular SSE: PCA, regression:",
      perpendicular_sse_pca, perpendicular_sse_regression)

grid = np.linspace(X[:, 0].min() - 0.5, X[:, 0].max() + 0.5, 200)

fig, ax = plt.subplots()
ax.scatter(X[:, 0], X[:, 1], label="data")
ax.plot(
    grid,
    xbar[1] + pca_slope * (grid - xbar[0]),
    label="PCA line",
)
ax.plot(
    grid,
    beta[0] + beta[1] * grid,
    label="regression line",
)
ax.scatter([xbar[0]], [xbar[1]], marker="x", s=80, label="mean")
ax.set_aspect("equal", adjustable="box")
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.legend()
plt.show()

Run and compare. The PCA line has slope (1), so its equation is [ y=x. ] The regression coefficients are [ b=0.4, \qquad m=0.8, ] so the regression line is [ \widehat y=0.4+0.8x. ]

The regression line has the smaller vertical sum of squares: [ 3.6<4. ] The PCA line has the smaller perpendicular sum of squares: [ 2<2.195\ldots. ]

Interpretation check. Each line is best for its own objective. The fact that both pass through the mean does not make their residuals or optimization problems the same.

Common mistake. In the data plot, regression residuals are vertical. Regression is an orthogonal projection in response-vector space, not perpendicular projection onto the displayed regression line.

Part 6. PCA from the SVD

Math reminder. If the centered data matrix has an SVD [ Z=U\Sigma V^T, ] then [ C=\frac1nZ^TZ = V\left(\frac1n\Sigma^T\Sigma\right)V^T ] and [ T=ZV=U\Sigma. ]

Thus the right singular vectors are principal directions, the squared singular values divided by (n) are covariance eigenvalues, and the columns of (U\Sigma) are score vectors.

Predict before running.

1. What singular values should correspond to covariance eigenvalues (9/2) and (1/2) when (n=4)?
2. How should the right singular directions compare with the covariance eigenvectors?

U_data, s_data, Vt_data = np.linalg.svd(Z, full_matrices=False)
V_data = Vt_data.T

covariance_eigenvalues_from_svd = s_data**2 / len(X)

T_from_svd = U_data @ np.diag(s_data)
T_direct = Z @ V_data

# Absolute values remove the irrelevant sign choices.
direction_alignment = np.abs(V_pca.T @ V_data)

Zhat1_svd = np.outer(T_from_svd[:, 0], V_data[:, 0])

pca_svd_scores_ok = np.allclose(T_from_svd, T_direct)
pca_svd_reconstruction_ok = np.allclose(Zhat1_svd, Zhat_diag)

fraction_first_svd = s_data[0]**2 / np.sum(s_data**2)

(
    s_data,
    covariance_eigenvalues_from_svd,
    direction_alignment,
    pca_svd_scores_ok,
    pca_svd_reconstruction_ok,
    fraction_first_svd,
    Zhat1_svd,
)

Run and compare. The singular values are approximately [ 3\sqrt2 \qquad\text{and}\qquad \sqrt2. ] Their squares divided by (4) are [ \frac92 \qquad\text{and}\qquad \frac12, ] the covariance eigenvalues.

The absolute direction-alignment matrix is the identity, up to roundoff. Thus the SVD and covariance calculations produce the same two axes, although their vector signs may differ.

Both score and reconstruction checks return True. The first-direction reconstruction is again [ \begin{bmatrix} 3/2&3/2\\ 3/2&3/2\\ -3/2&-3/2\\ -3/2&-3/2 \end{bmatrix}, ] and the captured fraction is (9/10).

Interpretation check. PCA from covariance and PCA from the SVD are two computations of the same principal directions, scores, and reconstructions.

Part 7. SVD, rank, and the four fundamental subspaces

Consider [ A= \begin{bmatrix} 1&-1\\ -2&2\\ 2&-2 \end{bmatrix}. ]

Its two columns are opposites, so its rank is (1).

Math reminder. For a full SVD [ A=U\Sigma V^T, ] the singular vectors associated with positive singular values give bases for the row and column spaces. The remaining right singular vectors give the null space, while the remaining left singular vectors give the left null space.

Predict before running.

1. What should the shapes of full (U), s, and (V^T) be for a (3\times2) matrix?
2. How many positive singular values should appear?
3. Which factor contains input-space directions? Which contains output-space directions?

A_map = np.array([[ 1.0, -1.0],
                  [-2.0,  2.0],
                  [ 2.0, -2.0]])

U_full, s_full, Vt_full = np.linalg.svd(
    A_map,
    full_matrices=True,
)

tolerance = 1e-10
rank_from_svd = int(np.sum(s_full > tolerance))

# Rows of Vt are transposed right singular vectors.
row_basis = Vt_full[:rank_from_svd, :]
null_basis = Vt_full[rank_from_svd:, :]

# Columns of U are left singular vectors.
column_basis = U_full[:, :rank_from_svd]
left_null_basis = U_full[:, rank_from_svd:]

svd_identity_ok = np.allclose(
    A_map @ Vt_full.T[:, 0],
    s_full[0] * U_full[:, 0],
)

null_ok = np.allclose(
    A_map @ null_basis.T,
    np.zeros((A_map.shape[0], null_basis.shape[0])),
)

left_null_ok = np.allclose(
    A_map.T @ left_null_basis,
    np.zeros((A_map.shape[1], left_null_basis.shape[1])),
)

print("shapes:", A_map.shape, U_full.shape, s_full.shape, Vt_full.shape)
print("singular values:", s_full)
print("rank:", rank_from_svd)
print("row-space basis (rows):\n", row_basis)
print("null-space basis (rows):\n", null_basis)
print("column-space basis (columns):\n", column_basis)
print("left-null-space basis (columns):\n", left_null_basis)
print("checks:", svd_identity_ok, null_ok, left_null_ok)

Run and compare. The shapes are [ A:\ (3,2), \qquad U:\ (3,3), \qquad s:\ (2,), \qquad V^T:\ (2,2). ] There is one positive singular value, so the numerical rank is (1).

The rows of row_basis and null_basis are vectors in the input space (\mathbb R^2). The columns of column_basis and left_null_basis are vectors in the output space (\mathbb R^3).

All three checks return True:

• the positive singular-vector identity holds;
• the right null-space basis is sent to zero by (A);
• the left-null-space basis is sent to zero by (A^T).

Interpretation check. A reduced SVD contains the factors needed to reconstruct (A), but a full (U) is useful here because its remaining columns display the full left null space.

Common mistake. The one-dimensional array s has length (\min(m,n)). The additional left-null directions are extra columns of full (U); they are not additional entries of s.

Part 8. Review and sanity checks

Before running the final cell, locate the earlier calculation corresponding to each statement.

1. One constraint gives parallel gradients.
2. Several constraints give a linear combination of constraint gradients.
3. The largest eigenvalue of (A^TA) is the squared maximum stretch.
4. Centering makes each feature column of (Z) have mean (0).
5. np.outer(scores, v) combines each scalar score with the direction (\mathbf v).
6. PCA wins the perpendicular-residual comparison, while regression wins the vertical-residual comparison.
7. The SVD of (Z) reproduces the PCA directions, scores, and reconstruction.
8. Positive and zero singular directions identify the four fundamental subspaces.

checks = {
    "one-constraint maximum": one_max_ok,
    "one-constraint minimum": one_min_ok,
    "several-constraint gradient span": several_ok,
    "sampled maximum stretch": np.isclose(
        sampled_stretches[j_max],
        singular_values_stretch[0],
    ),
    "sampled and eigenvector directions agree": np.isclose(
        sample_alignment,
        1.0,
    ),
    "PCA residuals are perpendicular": residuals_ok,
    "PCA Pythagorean identity": pythagoras_ok,
    "principal directions are orthonormal": orthogonal_ok,
    "principal coordinates diagonalize covariance": diagonal_ok,
    "score covariance matches": score_cov_ok,
    "regression wins vertical SSE": (
        vertical_sse_regression <= vertical_sse_pca + 1e-12
    ),
    "PCA wins perpendicular SSE": (
        perpendicular_sse_pca <= perpendicular_sse_regression + 1e-12
    ),
    "PCA scores agree with SVD": pca_svd_scores_ok,
    "PCA reconstruction agrees with SVD": pca_svd_reconstruction_ok,
    "singular-vector identity": svd_identity_ok,
    "right null-space check": null_ok,
    "left null-space check": left_null_ok,
}

assert all(checks.values())

checks

Run and compare. Every value should be True.

Exam-style checks.

1. Given a candidate and its constraint gradients, decide whether the code is checking one constraint or several constraints.
2. Given a covariance matrix, identify which eigenvector is PC1 and explain why unit length is required.
3. Given scores, Zhat, and R, distinguish coordinates, reconstructed vectors, and residuals.
4. Explain why changing the sign of a principal direction changes its scores but not its reconstructed data.
5. Given U, s, and Vt, identify input directions, output directions, stretch factors, rank, and the four fundamental subspaces.
6. Explain the identities [ \lambda_j=\frac{\sigma_j^2}{n}, \qquad T=ZV=U\Sigma. ]