Lab U6: SVD and Compression

Unit: Unit 6, Constraints and SVD
Role: Required
Textbook sections: maximum stretch and singular values; singular value decomposition; low-rank approximation and compression; applications and computation recap
Core path: SVD shapes, singular values, singular-vector identities, rank and null directions, redundant features, energy retained, rank-(k) reconstruction, compression on a small matrix, low-rank update shapes, and review code-reading checks
This lab is about reading short NumPy snippets as singular value decomposition. The goal is not to write long programs. The goal is to predict small outputs, identify shapes, translate code into mathematics, and explain what the computation says about the matrix map (\mathbf{x}\mapsto A\mathbf{x}).

Parts 0-9 are the core path. The final part contains review code-reading 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.3: Indexing and slicing
• Appendix B.4: Elementwise arithmetic versus linear algebra
• Appendix B.9: Numerical checks and roundoff
• Appendix B.10: NumPy linear algebra commands
• Appendix B.12: Small Python patterns used in the labs

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

Part 0. SVD shape habits

Math reminder. If (A) is (m\times n), then the reduced SVD in NumPy is read as [ A=U\Sigma V^T. ] The command

U, s, Vt = np.linalg.svd(A, full_matrices=False)

returns the left singular vectors in U, the singular values in the one-dimensional array s, and (V^T) in Vt. In the reduced SVD, U contains only the left singular-vector columns needed for the factorization, and s stores the diagonal entries of the reduced diagonal matrix (\Sigma).

Predict before running. If A is 3 x 2, what shapes should U, s, and Vt have when full_matrices=False?

import numpy as np

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

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

U, s, Vt = np.linalg.svd(A, full_matrices=False)

A.shape, U.shape, s.shape, Vt.shape, s

Run and compare. The output shapes are (3, 2), (3, 2), (2,), and (2, 2). The singular values are approximately [4.243, 0.].

Interpretation check. The array s stores only the diagonal entries of (\Sigma), not the full diagonal matrix. The zero singular value says that the matrix map forgets one input direction.

Common mistake. Vt is (V^T), not (V). To access right singular vectors as columns, use Vt.T.

Part 1. Diagonal stretch and singular values

Math reminder. Singular values are stretch factors. For a diagonal stretch, the stretch factors can often be read directly.

Predict before running. For [ A= \begin{bmatrix} 3&0\\ 0&1 \end{bmatrix}, ] which input direction should be stretched most?

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

e1 = np.array([1.0, 0.0])
e2 = np.array([0.0, 1.0])

A @ e1, A @ e2, np.linalg.norm(A @ e1), np.linalg.norm(A @ e2)

Run and compare. The outputs are [3, 0] and [0, 1]. The lengths are 3 and 1.

Interpretation check. The direction (\mathbf{e}_1) is stretched by (3). The direction (\mathbf{e}_2) is stretched by (1).

Review check. Which number is [ \max_{|\mathbf{x}|=1}|A\mathbf{x}|? ]

Part 2. (A^T A) and squared singular values

Math reminder. [ |A\mathbf{x}|^2=\mathbf{x}^TA^TA\mathbf{x}. ] Eigenvalues of (A^TA) are squared singular values.

Predict before running. For the matrix below, what should s**2 match?

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

B = A.T @ A
eigvals, eigvecs = np.linalg.eigh(B)

U, s, Vt = np.linalg.svd(A, full_matrices=False)

B, eigvals, s, s**2

Run and compare. The matrix (A^TA) is [ \begin{bmatrix} 5&4\\ 4&5 \end{bmatrix}. ] The eigenvalues are approximately 1 and 9, while s**2 is approximately 9 and 1.

Interpretation check. NumPy lists singular values from largest to smallest. The command np.linalg.eigh lists eigenvalues in increasing order. The same numbers appear, but in different orders.

Common mistake. Singular values are not eigenvalues of (A^TA). Their squares are eigenvalues of (A^TA).

Part 3. Checking (A\mathbf{v}_i=\sigma_i\mathbf{u}_i)

Math reminder. In an SVD, [ A\mathbf{v}_i=\sigma_i\mathbf{u}_i ] for each positive singular value.

Predict before running. Which variable below stores the first right singular vector?

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

U, s, Vt = np.linalg.svd(A, full_matrices=False)

i = 0
v = Vt.T[:, i]
lhs = A @ v
rhs = s[i] * U[:, i]

v, lhs, rhs, np.allclose(lhs, rhs)

Run and compare. The vector v is the first right singular vector. The arrays lhs and rhs agree, and np.allclose(lhs, rhs) returns True.

Interpretation check. The code checks [ A\mathbf{v}_1=\sigma_1\mathbf{u}_1. ] For this matrix, the stretch factor is (\sigma_1=3).

Common mistake. Vt[:, i] is not the (i)-th right singular vector. Use Vt.T[:, i].

Part 4. Rank and null directions from SVD

Math reminder. Zero singular values correspond to forgotten input directions. A right singular vector for a zero singular value lies in (\operatorname{null}(A)).

Predict before running. The two columns of the matrix below are opposites. What rank should the matrix have?

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

U, s, Vt = np.linalg.svd(A, full_matrices=False)

z = Vt[-1, :]

s, z, A @ z, np.linalg.matrix_rank(A)

Run and compare. The singular values are approximately [4.243, 0.]. The vector z is approximately [ \frac{1}{\sqrt2} \begin{bmatrix} 1\\ 1 \end{bmatrix}, ] up to sign, and A @ z is the zero vector.

Interpretation check. The matrix has rank (1). The input direction represented by z is forgotten.

Common mistake. Singular vectors are determined only up to sign. A vector and its negative represent the same singular direction.

Part 5. Redundant features revisited

Math reminder. In Unit 2, a redundant feature produced a nonzero null-space direction. SVD detects the same issue using a zero singular value.

Predict before running. The third column of (X) is the sum of the first two columns. What should happen to [ \begin{bmatrix} 1\\ 1\\ -1 \end{bmatrix} ] when multiplied by (X)?

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

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

U, s, Vt = np.linalg.svd(X, full_matrices=False)

X @ z, s, Vt[-1, :]

Run and compare. The output X @ z is zero. The singular values are approximately [8.478, 1.459, 0.]. The last row of Vt points in the same direction as [ \begin{bmatrix} 1\\ 1\\ -1 \end{bmatrix}, ] up to scaling and sign.

Interpretation check. The zero singular value detects the redundant-feature relation [ \text{column }1+\text{column }2-\text{column }3=\mathbf{0}. ]

Review check. If (X\mathbf{c}=\mathbf{y}), what is [ X(\mathbf{c}+5\mathbf{z})? ]

Part 6. Rank-(k) reconstruction

Math reminder. A rank-(k) reconstruction keeps the first (k) singular directions: [ A_k=U_k\Sigma_kV_k^T. ] In NumPy:

Ak = U[:, :k] @ np.diag(s[:k]) @ Vt[:k, :]

Predict before running. For the diagonal matrix below, what should a rank-2 reconstruction discard?

A = np.diag([3.0, 2.0, 0.5])

U, s, Vt = np.linalg.svd(A, full_matrices=False)

k = 2
Ak = U[:, :k] @ np.diag(s[:k]) @ Vt[:k, :]

s, Ak, np.linalg.matrix_rank(Ak)

Run and compare. The rank-2 reconstruction is [ \begin{bmatrix} 3&0&0\\ 0&2&0\\ 0&0&0 \end{bmatrix}. ]

Interpretation check. The reconstruction keeps the directions with singular values (3) and (2), and discards the direction with singular value (0.5).

Common mistake. The rank-(k) reconstruction is usually an approximation, not the original matrix.

Part 7. Energy retained

Math reminder. The energy retained by the first (k) singular values is [ \frac{\sigma_1^2+\cdots+\sigma_k^2} {\sigma_1^2+\cdots+\sigma_r^2}. ]

Predict before running. If the singular values are 10, 4, 1, 0.5, 0.1, should one singular value retain at least (90%) of the squared energy?

s = np.array([10.0, 4.0, 1.0, 0.5, 0.1])

energy = np.cumsum(s**2) / np.sum(s**2)

energy

Run and compare. The first singular value retains about 0.853 of the energy. The first two retain about 0.989.

Interpretation check. The smallest (k) retaining at least (90%) of the energy is (k=2).

Common mistake. Energy uses s**2, not s.

Part 8. Compression on a small matrix

Math reminder. A small grayscale image can be viewed as a matrix. Rank-(k) reconstruction keeps the strongest singular directions of that matrix. Here we use a tiny matrix instead of an image file.

Predict before running. The matrix below has one strong block and one weaker block. What should a rank-1 reconstruction keep?

M = np.array([[5.0, 5.0, 0.0, 0.0],
              [5.0, 5.0, 0.0, 0.0],
              [0.0, 0.0, 2.0, 2.0],
              [0.0, 0.0, 2.0, 2.0]])

U, s, Vt = np.linalg.svd(M, full_matrices=False)

M1 = U[:, :1] @ np.diag(s[:1]) @ Vt[:1, :]
M2 = U[:, :2] @ np.diag(s[:2]) @ Vt[:2, :]

s, np.cumsum(s**2) / np.sum(s**2), M1, M2

Run and compare. The singular values are approximately [10, 4, 0, 0]. The rank-1 reconstruction keeps the stronger (5)-block and discards the weaker (2)-block. The rank-2 reconstruction recovers the matrix exactly.

Interpretation check. Compression keeps dominant directions first. Smaller singular values represent weaker directions.

Review check. Why does the rank-2 reconstruction recover this particular matrix exactly?

Part 9. Low-rank update shapes

Math reminder. A low-rank update has the form [ W_{\text{new}}=W+BC, ] where [ B\in\mathbb R^{m\times r}, \qquad C\in\mathbb R^{r\times n}, ] and (r) is small.

Predict before running. If (m=1000), (n=800), and (r=4), what is the shape of (BC)? What is the largest possible rank of (BC)?

m, n, r = 1000, 800, 4

W_shape = (m, n)
B_shape = (m, r)
C_shape = (r, n)
update_shape = (B_shape[0], C_shape[1])

storage_update = m*r + r*n
storage_full = m*n

W_shape, B_shape, C_shape, update_shape, storage_update, storage_full

Run and compare. The update (BC) has shape (1000\times 800), so it can be added to (W). The matrices (B) and (C) store 7200 numbers, compared with 800000 numbers for a full (1000\times 800) matrix.

Interpretation check. Since the inner dimension is (4), [ \operatorname{rank}(BC)\leq 4. ]

Common mistake. A low-rank update does not mean that (W+BC) has low rank. It means the change (BC) has low rank.

Part 10. Review code-reading bank

Use these checks to practice reading short SVD snippets.

Check 1. Shape reading. In the next cell, identify the shapes of U[:, :k], np.diag(s[:k]), and Vt[:k, :]. Then compare them with the shapes produced by the wrong slice Vt[:, :k].

Check 2. Singular-vector identity. Which line checks [ A\mathbf{v}_1=\sigma_1\mathbf{u}_1? ]

Check 3. Debugging. Why would Vt[:, :k] be the wrong slice for rank-(k) reconstruction?

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

U, s, Vt = np.linalg.svd(A, full_matrices=False)

k = 1

good = U[:, :k] @ np.diag(s[:k]) @ Vt[:k, :]
identity_check = np.allclose(A @ Vt.T[:, 0], s[0] * U[:, 0])

good_shapes = (U[:, :k].shape, np.diag(s[:k]).shape, Vt[:k, :].shape)
bad_shapes = (U[:, :k].shape, np.diag(s[:k]).shape, Vt[:, :k].shape)

identity_check, good, good_shapes, bad_shapes

Run and compare. The identity check returns True. The rank-1 reconstruction is [ \begin{bmatrix} 3&0\\ 0&0 \end{bmatrix}. ] The tuple good_shapes gives the shapes used in the correct reconstruction: [ (2\times 1),\quad (1\times 1),\quad (1\times 2). ] The tuple bad_shapes shows the wrong final shape: [ (2\times 1),\quad (1\times 1),\quad (2\times 1). ]

Interpretation check. The slice Vt[:k, :] keeps the first (k) rows of (V^T), which are the transposed right singular vectors. The slice Vt[:, :k] keeps columns of (V^T), which is the wrong mathematical object and gives incompatible shapes in the reconstruction product.

Exam-style check. Given a short SVD computation, explain singular values, rank, retained energy, forgotten directions, and rank-(k) reconstruction without coding the SVD from scratch.