M3.6 · LESSON 6 OF 6

SVD and low rank

Multiply matrices, track dimensions and interpret transformations, inverses, eigenvectors and SVD.

THE IDEA

Every matrix is a rotation, a stretch, and a rotation

The singular value decomposition writes any matrix as A = U Σ Vᵀ: Vᵀ turns the input to line up with special directions, Σ stretches along each by a nonnegative singular value σ, and U turns the result into place. Equivalently, A is a sum of simple rank-one pieces: A = σ₁ u₁v₁ᵀ + σ₂ u₂v₂ᵀ + …, ordered from the biggest σ to the smallest. Each piece is one outer product, a column times a row.

If most of A's content sits in the first few pieces, you can keep only those. Keeping the top k gives the best possible rank-k approximation, and the error is exactly the square root of the sum of the squared discarded singular values.

Choose the rank yourself

Below is a tiny 10 by 12 grayscale picture: 120 numbers. Slide the rank k and watch the reconstruction, the error that is lost, and how many numbers you would store. The factors cost k·(10 + 12 + 1) numbers: k vectors of length 10, k of length 12, and k singular values.

original (10 by 12)

rank 1 approximation

what is lost

rank k = 1

kept singular values: 23.039

error = sqrt(sum of the discarded σ²) = 9.249 (measured: 9.249)

relative error = 37.254 percent of the image's norm

stored: 1 × (10 + 12 + 1) = 23 numbers, against 120 for the image

The factors are 80.833 percent smaller than the image.

Blue bars are the kept singular values. The picture is built from four simple pieces, so the later singular values are numerically zero.

The singular values here are 23.039, 6.766, 6.254, 0.802, 0, … Rank 1 already has an error of 37.254 percent; rank 3 3.232 percent; rank 4 is exact, for 92 stored numbers instead of 120.

Predict first. How many nonzero singular values can a rank-one matrix have?

When low rank helps, and when it does not

Low rank saves space only when k is small compared with the matrix: at rank 6 the factors need 138 numbers, more than the 120 you started with. The idea pays off for big matrices with fast-dropping singular values, such as compressing images or fine-tuning large models with a small low-rank change. Real data has noise, so singular values are rarely exactly zero: decide on a tolerance, or on how much of the total squared size you want to keep.

SVD with numpy

import numpy as np
U, s, Vt = np.linalg.svd(A, full_matrices=False)   # s sorted from large to small
A_k = (U[:, :k] * s[:k]) @ Vt[:k, :]               # best rank-k approximation
err = np.linalg.norm(A - A_k)                      # equals sqrt((s[k:] ** 2).sum())
rank = int((s > 1e-10 * s[0]).sum())               # count with a tolerance

Chapter context

A matrix represents a linear map once input and output bases are chosen. Multiplication composes maps. Tracking shapes is an algebraic requirement and also a practical way to detect implementation mistakes before running code.

TRY IT YOURSELF

Part 1. Do it by hand

The image is fixed. Your decisions are the rank k: how small an error you accept against how many numbers you store. Smallest k for 10 percent is not obvious until you try.

original (10 by 12)

rank 1 approximation

what is lost

rank k = 1

kept singular values: 23.039

error = sqrt(sum of the discarded σ²) = 9.249 (measured: 9.249)

relative error = 37.254 percent of the image's norm

stored: 1 × (10 + 12 + 1) = 23 numbers, against 120 for the image

The factors are 80.833 percent smaller than the image.

To do: Find the smallest rank k whose error is at most 10 percent of the image
To do: Find the smallest k that reproduces the image exactly (error below 0.001 percent)
To do: Find a k where storing the factors takes at least as many numbers as the image itself

Part 2. Do it in code

Write low_rank(A, k): the best rank-k approximation from np.linalg.svd, keeping the k largest singular values. Raise ValueError if k is less than 1 or larger than the smaller dimension of A.
Write error_after(A, k): the Frobenius norm of A minus low_rank(A, k). It should equal the square root of the sum of the squared discarded singular values.
Write numerical_rank(A, tol=1e-10): the number of singular values larger than tol times the largest one. Count by singular values, not by exact zeros.

import numpy as np

def low_rank(A, k):
    # keep the k largest singular values
    pass

def error_after(A, k):
    # Frobenius norm of A - low_rank(A, k)
    pass

def numerical_rank(A, tol=1e-10):
    # how many singular values are above tol times the largest
    pass

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

Ctrl/⌘ + Enter runs it

Stuck, or something not working? You can come back to it any time.

CHECK YOUR UNDERSTANDING

How many nonzero singular values can a rank-one matrix have?

ANone; a rank-one matrix has all singular values equal to zero.BTwo; one for each of the two factors in the outer product.CExactly one.

Chapter reference: math, worked example & code

This section belongs to the whole chapter, so it is the same on every lesson in it and may cover a different lesson's topic. The code for the experiment you just ran is inside the experiment itself.

Math and notation

C_{ij}=\sum_k A_{ik}B_{kj},\quad (AB)^T=B^TA^T.

\det\begin{bmatrix}a&b\\c&d\end{bmatrix}=ad-bc,\quad Av=\lambda v.

A=U\Sigma V^T.

U and V contain orthonormal singular directions; diagonal entries of Sigma are singular values. Eigenvalues and singular values are different concepts.

For an invertible 2x2 matrix, the inverse is

\begin{bmatrix}a&b\\c&d\end{bmatrix}^{-1}=\frac1{ad-bc}\begin{bmatrix}d&-b\\-c&a\end{bmatrix}.

This formula requires a nonzero determinant. A near-zero determinant can signal numerical sensitivity, so exact invertibility alone does not guarantee a well-conditioned computation.

Worked example

Multiply A=[[1,2],[3,4]] by B=[[2,0],[1,2]]. Row-column products give C=[[4,4],[10,8]]. The determinant of A is 1×4−2×3=−2, so A is invertible. In general AB is not equal to BA; their shapes may even differ.

Reference code

The following example uses the Python standard library and does not call a remote service.

A, B = [[1,2],[3,4]], [[2,0],[1,2]]
C = [[sum(A[i][k]*B[k][j] for k in range(2)) for j in range(2)] for i in range(2)]
assert C == [[4,4],[10,8]]
assert A[0][0]*A[1][1] - A[0][1]*A[1][0] == -2

Python examples are provided for study and local execution; this page does not execute Python.

Common mistakes

AB usually differs from BA. A transpose is not generally an inverse. A negative eigenvalue is possible, while singular values are nonnegative.