SVD: rotate, stretch, rotate

In Lessons 11 and 12, eigenvectors were a very handy tool. Switch to an eigenvector basis and the matrix becomes "just stretching", and powers become easy. But there was a catch. A rotation matrix had no real eigenvectors, and a shear matrix had only one. On top of that, eigenvectors are only defined for square matrices. You cannot use them on a data table like 300 users × 1000 movies.

This lesson's singular value decomposition (SVD) clears up all of that. Any matrix, square or not, can always be split into three steps: rotate → stretch along the axes → rotate again. Having no exceptions is what made SVD the most widely used decomposition in linear algebra.

It shows up everywhere, too. Compressing a photo into fewer numbers, removing noise, finding the "taste axes" of a recommender system, finding the null space (which Lesson 9 only used as a tool), and next lesson's PCA all stand on SVD.

14-1. A circle always becomes an ellipse (5 min)

In Lesson 5, a matrix was a rule that moves the whole grid. This time, instead of the grid, we move the unit circle (every point at distance 1 from the origin).

Whatever 2×2 matrix you use, the unit circle always becomes an ellipse. Sometimes it stays a circle (a rotation), sometimes it gets squashed all the way flat into a line segment (det = 0), but it never turns into a lumpy potato shape. That is because the transformation keeps grid lines parallel and evenly spaced.

An ellipse has a longest radius and a shortest radius, and the two directions are perpendicular. SVD is the story of finding those two radii.

• The length of the ellipse's long radius = the first singular value σ₁
• The length of the ellipse's short radius = the second singular value σ₂
• The points on the circle that go to those two radii (perpendicular unit vectors) = the right singular vectors v₁, v₂
• The directions of the ellipse's two radii (also perpendicular) = the left singular vectors u₁, u₂

The key point is that two perpendicular arrows on the circle are still perpendicular after being moved. An ordinary matrix does not keep right angles (move î and ĵ with a shear matrix and they end up slanted). But on the circle there is always one special pair whose right angle is kept. That pair is v₁, v₂.

14-2. Rotate, stretch, rotate (10 min)

In the lab below, the blue (indigo) shape starts as the unit circle. Press the buttons on the right in order.

• Start (unit circle): the circle and two yellow arrows. The dark yellow one is v₁ and the pale yellow one is v₂.
• ① Vᵀ: rotate: rotates the circle so that v₁ lies on the x-axis and v₂ on the y-axis. A rotated circle is still a circle, so the shape does not change.
• ② Σ: stretch along the axes: stretches by σ₁ along the x-axis and by σ₂ along the y-axis. Now the circle is an ellipse lined up with the axes.
• ③ U: rotate again = A: rotates the whole ellipse into its final place. The result is exactly what you get by moving the circle with the matrix A.

The yellow arrows follow where the right singular vectors go. In the last step, check that the tips of the yellow arrows land exactly on the ends of the ellipse's long and short radii.

Things worth trying in the lab. The handles for the tips of A's columns only appear at step ③ (= A).

• Put column 1 at (0, 1) and column 2 at (−1, 0) and you have a 90° rotation. The circle stays a circle, so σ₁ = σ₂ = 1. That means the matrix does no stretching at all.
• Put column 1 at (1, 0) and column 2 at (1, 1) and you have the shear from Lesson 6. In Lesson 11 it gave us trouble by having only one eigenvector, but SVD handles it without a fuss: σ₁ ≈ 1.618, σ₂ ≈ 0.618. The golden ratio from Lesson 12 shows up again here.
• Swap columns 1 and 2 to make det negative, and at step ③ the ellipse squashes flat once and comes out flipped. In that case U is not a pure rotation but a rotation with a flip mixed in (Lesson 7: a negative det means the orientation flips).

Reading the book's notation

Writing the three steps in one line gives the singular value decomposition. Matrix products are read from the right (Lesson 6), so $V^{\mathsf T}$ comes first and $U$ comes last.

The rightmost $V^{\mathsf T}$ is ① rotate, the $\Sigma$ in the middle is ② stretch, and the $U$ on the left is ③ rotate again. The columns of U are $\mathbf{u}_1, \mathbf{u}_2$, and the rows of $V^{\mathsf T}$ are $\mathbf{v}_1, \mathbf{v}_2$ laid on their sides.

Follow a single vector through the three steps and you get the most important formula of SVD. $\mathbf{v}_1$ becomes $(1, 0)$ on the x-axis in ①, then $(\sigma_1, 0)$ in ②, then $\sigma_1\mathbf{u}_1$ in ③.

Compare it with $A\mathbf{v} = \lambda\mathbf{v}$ from Lesson 11. An eigenvector had to come back "in the same direction", but a singular vector is allowed to go somewhere else ($\mathbf{u}$). In exchange, the two starting directions must be perpendicular, and so must the two arriving directions. Changing the condition this way gave every matrix an answer.

We called U and V "rotations", but their precise name is orthogonal matrices. Like the Q from Lesson 13's QR decomposition, their columns are perpendicular and have length 1, and they are either rotations or rotations with a flip mixed in. Either way, they do not change lengths or angles.

14-3. Singular values, eigenvalues, the determinant, and rank (6 min)

Singular values are tightly connected to the ideas you learned earlier.

The link to eigenvalues. Write $A^{\mathsf T}A$ out using the SVD, and the $U^{\mathsf T}U$ in the middle becomes $I$ and disappears.

This is the same shape as Lesson 12's diagonalization $PDP^{-1}$ ($V^{\mathsf T}$ is exactly $V^{-1}$). So the eigenvectors of $A^{\mathsf T}A$ are the $\mathbf{v}_i$, and its eigenvalues are $\sigma_i^2$. The lab computes σ and V this way too. $A^{\mathsf T}A$ is always symmetric, so as you saw in Lesson 11 it has perpendicular eigenvectors, and its eigenvalues are 0 or more. That is why every matrix has an SVD.

The link to the determinant. In Lesson 7, det was the area scale factor. Rotations do not change area; only the stretch scales area, by $\sigma_1\sigma_2$. So

The sign (whether it flipped) is carried by U and V.

The link to rank. The rank from Lesson 9 is "the number of nonzero singular values". If $\sigma_2 = 0$, step ② squashes the whole y direction away, and the circle becomes a line segment. Then $A\mathbf{v}_2 = 0\cdot\mathbf{u}_2 = \mathbf{0}$, so $\mathbf{v}_2$ is exactly the direction of the null space. That is why we found the null space with Vt[-1] in Lesson 9.

14-4. SVD in NumPy (12 min)

Let's follow the code in the "Check it with NumPy" panel below the lab from the beginning. We start with the lab's initial matrix A.

>>> import numpy as np
>>> A = np.array([[2, 1], [0.5, 1.5]])
>>> U, S, Vt = np.linalg.svd(A)
>>> S
array([2.55833637, 0.97719754])

np.linalg.svd returns three results. Line up three names with commas, as in U, S, Vt = …, and you receive them one at a time (the same syntax as vals, vecs = np.linalg.eig(A) in Lesson 11). Watch the names.

• S is not the matrix Σ but a 1-D array holding only the singular values, largest first.
• Vt is not V but Vᵀ, already transposed. So the right singular vector v₁ comes out as Vt[0], a row. The left singular vectors of U come out as U[:, 0], a column.

>>> U
array([[-0.85065081, -0.52573111],
       [-0.52573111,  0.85065081]])
>>> Vt
array([[-0.76775173, -0.64074744],
       [-0.64074744,  0.76775173]])

The NumPy panel below the lab rounds each line so the decimals do not run on too long: np.round(U, 6), np.round(S, 6), np.round(Vt, 6), np.round(U @ np.diag(S) @ Vt, 6) (that is why the panel's S reads array([2.558336, 0.977198])). np.round(array, decimals) gives a new array with every entry rounded to that many decimal places. It does not change the values; it is just for display.

>>> np.round(U, 3)
array([[-0.851, -0.526],
       [-0.526,  0.851]])

Multiply the three pieces back together and A should come back. Since S is a 1-D array, first build the diagonal matrix Σ with np.diag(S). Given a 1-D array, np.diag builds a square matrix with those numbers on its diagonal.

>>> np.diag(S)
array([[2.55833637, 0.        ],
       [0.        , 0.97719754]])
>>> U @ np.diag(S) @ Vt
array([[2. , 1. ],
       [0.5, 1.5]])
>>> np.allclose(U @ np.diag(S) @ Vt, A)
True

It came out clean this time, but floating-point arithmetic often leaves crumbs like 1.9999999999999998. That is why we compare with np.allclose, which we have used since Lesson 3.

Compare A @ Vt[0] with S[0] * U[:, 0] and you can check the formula Av₁ = σ₁u₁ from 14-2 directly.

>>> A @ Vt[0]
array([-2.1762509 , -1.34499702])
>>> S[0] * U[:, 0]
array([-2.1762509 , -1.34499702])
>>> np.allclose(A @ Vt[1], S[1] * U[:, 1])
True

Let's also check that U and V are orthogonal matrices, exactly the way we checked Q in Lesson 13.

>>> np.allclose(U.T @ U, np.eye(2))
True
>>> np.allclose(Vt @ Vt.T, np.eye(2))
True

Why the signs differ from the lab

The lab's U is roughly [[0.85, -0.53], [0.53, 0.85]], while NumPy's U is [[-0.85, -0.53], [-0.53, 0.85]]. The first column has the opposite sign. The first row of Vᵀ also has the opposite sign. Neither one is wrong.

If you flip v₁ and u₁ together in Av₁ = σ₁u₁, you get A(−v₁) = σ₁(−u₁), and the formula still holds. In other words, flipping a u and its partner v together gives the same decomposition. The lab picks signs so that V is a pure rotation, while NumPy (LAPACK underneath) hands back whatever signs its computation produced. Let's flip the pair together to match the lab.

>>> U2 = U.copy()
>>> Vt2 = Vt.copy()
>>> U2[:, 0] *= -1
>>> Vt2[0] *= -1
>>> np.round(U2, 2)
array([[ 0.85, -0.53],
       [ 0.53,  0.85]])
>>> np.round(Vt2, 2)
array([[ 0.77,  0.64],
       [-0.64,  0.77]])
>>> np.allclose(U2 @ np.diag(S) @ Vt2, A)
True

.copy() makes a fresh copy of the array. If you write U2 = U without copying, you only give the same array a second name, and changing U2 changes U too. U2[:, 0] *= -1 means "multiply all of column 0 by −1 and put it back in place".

Flip only one side and the decomposition breaks. When you compare two SVD results, compare the singular values S, or compare the vectors ignoring sign (np.abs).

The links to eigenvalues and the determinant, in numbers

>>> vals, vecs = np.linalg.eigh(A.T @ A)
>>> vals
array([0.95491503, 6.54508497])
>>> S ** 2
array([6.54508497, 0.95491503])
>>> np.linalg.det(A)
np.float64(2.5)
>>> S.prod()
np.float64(2.5000000000000004)

The eigenvalues of A.T @ A equal σ². The only difference is that eigh lists them smallest first and svd largest first. S.prod() multiplies all the numbers in the array together, and it differs from det 2.5 only by crumbs in the last digit.

Let's also try the shear matrix and the flat matrix from Lesson 9. When you only need the singular values, pass compute_uv=False and you get back just S, without computing U and V.

>>> np.linalg.svd(np.array([[1, 1], [0, 1]]), compute_uv=False)
array([1.61803399, 0.61803399])
>>> F = np.array([[2, 1], [1, 0.5]])
>>> U, S, Vt = np.linalg.svd(F)
>>> np.round(S, 6)
array([2.5, 0. ])
>>> np.linalg.matrix_rank(F)
np.int64(1)
>>> np.allclose(F @ Vt[-1], 0)
True

The flat matrix F has a second singular value of 0, rank 1, and sends Vt[-1] (v₂) to the zero vector. Depending on your setup, the second singular value may come out as a tiny number like 1e-17, so we rounded it for display. np.linalg.matrix_rank also runs an SVD inside and counts tiny singular values as 0.

Non-square matrices

SVD works for non-square matrices too. Feed it a 3×2 matrix and the three pieces have different shapes.

>>> B = np.array([[1, 0], [1, 1], [0, 1]])
>>> U, S, Vt = np.linalg.svd(B)
>>> U.shape, S.shape, Vt.shape
((3, 3), (2,), (2, 2))
>>> U @ np.diag(S) @ Vt
Traceback (most recent call last):
  ...
ValueError: matmul: Input operand 1 has a mismatch in its core dimension 0, with gufunc signature (n?,k),(k,m?)->(n?,m?) (size 2 is different from 3)

By default (full_matrices=True), U comes out as a full 3×3 orthogonal matrix. But there are only 2 singular values, so U's third column has nothing to pair with, and its shape does not fit np.diag(S) (2×2), hence the error. To get only the columns that are actually used, pass full_matrices=False.

>>> U, S, Vt = np.linalg.svd(B, full_matrices=False)
>>> U.shape, S.shape, Vt.shape
((3, 2), (2,), (2, 2))
>>> np.allclose(U @ np.diag(S) @ Vt, B)
True

For matrices with very many rows, like data tables, full_matrices=False also uses far less memory. Next lesson's PCA uses it this way too.

14-5. Low-rank approximation: keep only the large singular values (10 min)

Write $A = U\Sigma V^{\mathsf T}$ out by columns and rows, and the matrix becomes a sum of rank-1 pieces.

$\mathbf{u}\mathbf{v}^{\mathsf T}$ is the matrix you get by multiplying a column vector by a row vector (the outer product). Every row is a multiple of $\mathbf{v}^{\mathsf T}$, so its rank is 1. Each piece carries a singular value like a weight, and the singular values are lined up largest first. So if you keep only the first few pieces and throw the rest away, you can write a matrix very close to the original with far fewer numbers. This is called the rank-k approximation.

Just as Lesson 13's least-squares line was the line "no other line can beat on residuals", $A_k$ is the matrix closest to A among all matrices of rank k (the Eckart–Young theorem). The pieces you throw away are exactly the error, so you can read the size of the error straight off the discarded singular values.

Here $\lVert\cdot\rVert$ means: square every entry of the matrix, add them up, and take the square root. It is the norm from Lesson 3 stretched to matrices.

First, split the 2×2 matrix A into its two pieces. np.outer(u, v) is the outer product u vᵀ.

>>> U, S, Vt = np.linalg.svd(A)
>>> A1 = S[0] * np.outer(U[:, 0], Vt[0])
>>> A2 = S[1] * np.outer(U[:, 1], Vt[1])
>>> np.round(A1, 3)
array([[1.671, 1.394],
       [1.033, 0.862]])
>>> np.linalg.matrix_rank(A1)
np.int64(1)
>>> np.allclose(A1 + A2, A)
True

The first piece A1 alone already captures the broad shape of A = [[2, 1], [0.5, 1.5]] fairly well. In lab terms, A1 is the transformation that squashes the ellipse's short radius to 0, leaving a line segment.

Compressing a small grayscale photo

A grayscale photo is a matrix with one brightness per pixel (0 = black, 1 = white). Let's make a 32×32 photo ourselves: a bright square in the middle, a gray band across it, and small random values added, like camera noise.

>>> rng = np.random.default_rng(0)
>>> img = np.zeros((32, 32))
>>> img[8:24, 8:24] = 1.0
>>> img[14:18, 2:30] = 0.5
>>> clean = img.copy()
>>> img = img + 0.05 * rng.standard_normal((32, 32))
>>> img.shape
(32, 32)

img[8:24, 8:24] = 1.0 fills the square of rows 8–23 and columns 8–23 with 1.0 in one go (Lesson 1's slicing, used on two axes). rng.standard_normal((32, 32)) builds a 32×32 matrix of random numbers with mean 0 and standard deviation 1. Because we fixed the seed at 0, you get the same photo no matter how many times you run it. We kept a separate copy of the clean photo, before noise, in clean.

Now run the SVD and look at the singular values.

>>> U, S, Vt = np.linalg.svd(img)
>>> np.round(S[:6], 2)
array([14.45,  3.35,  0.49,  0.47,  0.45,  0.44])

Of the 32 singular values, only the first two are large; from the third on they are all similarly small. That means two patterns, the square and the band, are almost the whole photo, and the rest is noise. Let's define a function that builds the rank-k approximation.

>>> def approx(k):
...     return U[:, :k] @ np.diag(S[:k]) @ Vt[:k]
>>> approx(2).shape
(32, 32)
>>> np.linalg.matrix_rank(approx(2))
np.int64(2)

U[:, :k] is the first k columns of U, and Vt[:k] is the first k rows of Vᵀ. It is (32×k) @ (k×k) @ (k×32), so the result is a 32×32 photo again, but of rank k.

Let's vary k and compare the error against the original photo. The error is "the norm of the difference ÷ the norm of the photo", so 0.1 means roughly 10% off.

>>> for k in [1, 2, 3, 5, 10, 32]:
...     err = np.linalg.norm(img - approx(k)) / np.linalg.norm(img)
...     print(k, round(float(err), 3))
1 0.245
2 0.098
3 0.092
5 0.081
10 0.058
32 0.0

np.linalg.norm is the function that measured vector lengths in Lesson 3; give it a matrix and it returns the "square root of the sum of squares of all entries" described above. The error drops sharply going from k = 1 to 2, and only a little after that, just as the list of singular values suggested. Let's also check the formula saying the error is set by the discarded singular values.

>>> np.isclose(np.linalg.norm(img - approx(2)), np.sqrt(np.sum(S[2:] ** 2)))
np.True_

np.True_ is the true value NumPy hands back. Use it exactly like Python's True.

Count the numbers you have to store and the compression shows. The original photo is 32 × 32 = 1024 numbers. The rank-2 approximation can be rebuilt from just 2 columns of U (32 each), 2 rows of Vᵀ (32 each), and 2 singular values.

>>> 32 * 32
1024
>>> 2 * (32 + 32 + 1)
130

That is one-eighth of the numbers, with an error of about 10%. Real photos have many more patterns, so k is set in the tens, but the principle is the same.

Throwing away makes it cleaner

This time, compare the approximations with clean, the photo before noise was added.

>>> round(float(np.linalg.norm(img - clean)), 3)
1.557
>>> for k in [1, 2, 3, 5, 10]:
...     print(k, round(float(np.linalg.norm(approx(k) - clean)), 3))
1 3.347
2 0.542
3 0.73
5 0.975
10 1.298

The noisy photo img differs from the clean one by 1.557. Yet the rank-2 approximation is off by only 0.542: it is closer to the clean photo than the noisy original is. The pieces with small singular values mostly hold noise, so throwing them away sweeps the noise out with them. Raise k further and the noise comes back, moving you farther away again. Pull out a single row and it is even clearer.

>>> np.round(img[16, :8], 2)
array([-0.03, -0.01,  0.56,  0.45,  0.6 ,  0.59,  0.41,  0.49])
>>> np.round(approx(2)[16, :8], 2)
array([0.02, 0.  , 0.5 , 0.49, 0.52, 0.5 , 0.47, 0.47])

These are the first 8 pixels of row 16, where the band runs. The clean values are 0, 0, 0.5, 0.5, …, and the rank-2 approximation sits much more evenly around 0.5. Approximating a users × movies ratings table with a low rank in a recommender system is the same idea. A few large singular values become taste axes like "how much you like action", and the rest is treated as random wobble.

14-6. Common mistakes (3 min)

Transposing Vt one more time. The third result of np.linalg.svd is already Vᵀ. If you receive it as V and then write U @ np.diag(S) @ V.T, you get the wrong matrix. Get in the habit of naming it Vt from the start. Right singular vectors are Vt[i] (rows); left singular vectors are U[:, i] (columns).

Multiplying S as if it were a matrix. S is a 1-D array, so U @ S @ Vt is not multiplying by Σ. When you put it in a matrix product, turn it into np.diag(S). (U * S happens to multiply each column of U by its singular value through Lesson 2's broadcasting, giving the same values as U @ np.diag(S), but that is easy to mix up at first, so np.diag is the safer choice.)

Thinking a different sign means a wrong answer. The U and V from the lab, from NumPy and from other libraries can differ in sign. If a u and its partner v are flipped together, it is the same decomposition.

Exercises (5 min)

1. Predict the singular values of D = np.array([[3, 0], [0, -2]]), then check. Where did the −2 on the diagonal go?
2. In the lab, if you put column 1 at (2, 0) and column 2 at (0, 0.5), what are σ₁ and σ₂? Does the circle rotate in steps ① and ③?
3. For the photo in 14-5, how many numbers do you need to store for the rank-1 approximation alone? And what was its error against the original, in percent?

Solutions

>>> D = np.array([[3, 0], [0, -2]])
>>> U, S, Vt = np.linalg.svd(D)
>>> S
array([3., 2.])
>>> Vt[1, 1]
np.float64(-1.0)
>>> np.allclose(U @ np.diag(S) @ Vt, D)
True
>>> 1 * (32 + 32 + 1)
65

Question 1: singular values are stretch lengths, so they are never negative. They are 3 and 2, and the minus sign of −2 (the flip of the y-axis) was taken by Vᵀ. Indeed, Vt[1, 1] is −1.

Question 2: the matrix already only stretches along the axes, so σ₁ = 2 and σ₂ = 0.5, and V and U are matrices that do not rotate (the identity), so the circle does not move in ① and ③. Singular values must be lined up largest first, so if you flip it around and put column 1 at (0.5, 0) and column 2 at (0, 2), you can watch the circle rotate by 90° in ① and in ③.

Question 3: 32 + 32 + 1 = 65 numbers, and in the table in 14-5 the error was 0.245, about 25%.

In the next lesson you will use SVD and eigenvalues on a cloud of data points to find the direction of widest spread, which is principal component analysis (PCA), and look back over the road we have traveled in 15 lessons.