Singular value decomposition (SVD)

Every matrix — whatever its shape — can be written as three simpler operations in sequence:

$A = U\Sigma V^\top$

A matrix multiplication therefore always is: turn, stretch, turn. Nothing else.

The singular values in $\Sigma$ are sorted in order of size and say how much the matrix stretches in each direction. If the fifth value is small it means that the matrix barely does anything in the fifth direction — and it can then be thrown away almost for free.

That is the whole idea of low-rank approximation: keep the $k$ largest singular values, discard the rest.

For $A \in \mathbb{R}^{m\times n}$ of rank $r$:

$A = \sum_{i=1}^{r} \sigma_i\, u_i v_i^\top, \qquad \sigma_1 \geq \sigma_2 \geq \dots \geq \sigma_r > 0$

The Eckart–Young theorem: the best rank-$k$ approximation of $A$ (in the Frobenius and the spectral norm) is obtained simply by truncating the sum:

$A_k = \sum_{i=1}^{k} \sigma_i u_i v_i^\top, \qquad \|A - A_k\|_F^2 = \sum_{i>k}\sigma_i^2$

That is a remarkable result: the optimal approximation requires no search, it is read off directly.

The compression: $A$ has $mn$ numbers, $A_k$ has $k(m + n + 1)$. For a 1000×1000 matrix with $k = 50$ that is 100 050 against 1 000 000 — a tenth.

Four uses:

The connection to eigenvalues: $\sigma_i^2$ are the eigenvalues of $A^\top A$, and $v_i$ its eigenvectors. But SVD exists for every matrix, including non-square and singular ones — unlike the eigendecomposition. That is why it is the workhorse.

Computing it costs $O(mn\min(m,n))$ for the full SVD. If you only need the $k$ largest there are randomised methods that are dramatically faster — sklearn.utils.extmath.randomized_svd or scipy.sparse.linalg.svds.

import numpy as np

rng = np.random.default_rng(0)

# A matrix that REALLY has rank 3, plus a little noise
A = rng.normal(size=(200, 5)) @ rng.normal(size=(5, 150))
A = A[:, :3] @ rng.normal(size=(3, 150)) + 0.1 * rng.normal(size=(200, 150))

U, s, Vt = np.linalg.svd(A, full_matrices=False)
print(np.round(s[:8], 2))
# [419.4  406.63 242.5    2.61   2.5    2.47   2.44   2.42]
#   ↑ three large ones, then a jump down to the noise level

# The explained variance per component
share = s**2 / (s**2).sum()
print(np.round(np.cumsum(share)[:5], 4))     # [0.4394 0.8524 0.9993 0.9993 0.9993]

def truncate(U, s, Vt, k):
    return U[:, :k] * s[:k] @ Vt[:k]

for k in (1, 3, 10, 50):
    Ak = truncate(U, s, Vt, k)
    error = np.linalg.norm(A - Ak) / np.linalg.norm(A)
    stored = k * (A.shape[0] + A.shape[1] + 1)
    print(f"k={k:>3}  relative error {error:.4f}  storage {stored / A.size:.1%} of the original")
# k=  1  relative error 0.7487  storage 1.2% of the original
# k=  3  relative error 0.0268  storage 3.5% of the original
# k= 10  relative error 0.0248  storage 11.7% of the original

# Eckart–Young: the error is exactly the sum of the discarded squared singular values
k = 3
print(round(float(np.linalg.norm(A - truncate(U, s, Vt, k))**2), 4),
      round(float((s[k:]**2).sum()), 4))            # the same number

# The LoRA idea: a rank-r update of a large layer
d, r = 4096, 8
print(f"a full layer {d*d:,} parameters, LoRA rank {r}: {2*d*r:,} "
      f"({2*d*r/(d*d):.2%})")
# a full layer 16,777,216 parameters, LoRA rank 8: 65,536 (0.39%)