Eigenvectors, SVD and low rank

Problem

A weight matrix with 4096 × 4096 entries is 16.8 million numbers, but that count says little about what the matrix does. Some questions need more. Which input directions does it amplify, and which does it ignore? Can it be stored in far fewer numbers without changing its behaviour much? Will solving with it amplify a rounding error into a wrong answer? A matrix product is easy to compute but hard to read.

Decompositions answer these questions. They rewrite a matrix as a product of simple pieces: rotations, which keep lengths, and diagonal scalings, which you can read off. The same pieces explain why PCA works, why fine-tuning with LoRA can update a 4096 × 4096 layer with 65,536 parameters, and why some layers train badly.

Intuition

Push every point of the unit circle through a 2 × 2 matrix and you get an ellipse. That is true of every matrix. The ellipse's longest axis is the direction the matrix stretches most, and its length is the largest singular value σ₁. Its shortest axis has length σ₂. The input directions that land on those axes are the right singular vectors v₁, v₂, and the directions of the axes themselves are the left singular vectors u₁, u₂. So the matrix rotates v₁, v₂ onto the coordinate axes, stretches by σ₁, σ₂, and rotates the result onto u₁, u₂: A = U Σ Vᵀ.

A = [[3, 0], [4, 5]] sends v₁ to σ₁u₁ and v₂ to σ₂u₂. Every matrix has this shape: a rotation (or reflection), a stretch along the axes, and another rotation. Its largest stretch is 6.708, its smallest is 2.236, and their ratio is 3.0.

An eigenvector is the special case where no rotation happens: the matrix only stretches it, Av = λv. Symmetric matrices, such as covariances and Hessians, have a full set of eigenvectors at right angles to each other, and for them the two decompositions agree. Most other matrices, including nearly every weight matrix, rotate every direction somewhat. The SVD is the tool that works for all of them.

Mechanics

Eigen-decomposition

A v = λ v  ·  symmetric A = Q Λ Qᵀ,   QᵀQ = I

For a symmetric matrix S, the eigenvalues λᵢ are real, the eigenvectors can be chosen orthonormal, and S = Σᵢ λᵢ qᵢqᵢᵀ is a sum of rank-one pieces. For S = [[2, 1], [1, 2]] the eigenvalues are 1 and 3, with eigenvectors [0.707, -0.707] and [0.707, 0.707]. The matrix stretches the diagonal direction by 3 and the anti-diagonal by 1. torch.linalg.eigh handles symmetric (Hermitian) input, reads only one triangle (UPLO='L' by default), and returns eigenvalues in ascending order (eigh, eigvalsh). torch.linalg.eig takes any square matrix and always returns complex output, because a rotation has no real eigenvectors (eig).

Eigenvalues tell you how a matrix behaves when applied over and over. Aᵏv = λᵏv, so directions with |λ| > 1 blow up and those with |λ| < 1 die out. That is the heart of exploding and vanishing gradients in recurrent networks. Repeated multiplication also finds the top eigenvector: power iteration multiplies a vector by A and renormalizes it, and the component along the top eigenvector wins by a factor of λ₁/λ₂ each step.

The singular value decomposition

A = U Σ Vᵀ = Σᵢ σᵢ uᵢ vᵢᵀ,   A: [m, n], σ₁ ≥ σ₂ ≥ … ≥ 0

U's columns are orthonormal in the output space and V's in the input space, and Σ is diagonal. The singular values are the square roots of the eigenvalues of AᵀA, and the vᵢ are its eigenvectors. That is a good way to think about the SVD but a bad way to compute it, because forming AᵀA squares every singular value and loses small ones to rounding (see Tradeoffs). torch.linalg.svd returns (U, S, Vh) with Vh = Vᵀ and S in descending order. By default full_matrices=True pads U or Vh to square; pass False for the thin version you almost always want (svd). svdvals returns only S, and its gradients are always well defined (svdvals). On CUDA the default driver is gesvdj, a Jacobi method, which falls back to the QR-based gesvd if it fails.

Rank and low-rank approximation

The rank is the number of non-zero singular values: how many independent directions the matrix can output. In floating point nothing is exactly zero, so torch.linalg.matrix_rank counts singular values above max(atol, rtol · σ₁). When neither tolerance is given, rtol = ε · max(m, n), with ε the machine epsilon of the dtype (matrix_rank).

A_k = Σ_i≤k σᵢ uᵢ vᵢᵀ  ·  ‖A − A_k‖²_F = Σ_i>k σᵢ² = min_rank(B)≤k ‖A − B‖²_F

Eckart–Young says that cutting the SVD after k terms gives the best rank-k approximation, in both the Frobenius and the spectral norm. The error is exactly the singular values you dropped. A rank-k matrix stores k(m + n) numbers instead of mn, so when the singular values fall off fast you can compress a lot and lose little. For very large matrices, torch.svd_lowrank estimates only the top q singular triplets with a randomized method (Halko, Martinsson and Tropp). It multiplies A by a random [n, q] matrix, orthonormalizes the result with QR, refines it with niter=2 rounds of subspace iteration, and runs an exact SVD on the small q × n matrix QᵀA (get_approximate_basis, _svd_lowrank). Unlike linalg.svd, it returns V, not Vᵀ (svd_lowrank).

Two singular values carry almost everything. The last two segments are too thin to label: 0.6% and 0.6%. Each rank-k error equals the sum of the σᵢ² it drops.

The condition number

κ(A) = σ_max / σ_min  ·  ‖Δx‖/‖x‖ ≤ κ(A) · ‖Δb‖/‖b‖   for Ax = b

The condition number is the ratio of the ellipse's longest axis to its shortest. If κ is 10ᵏ, solving with A can lose about k digits: a relative change in b can grow by up to a factor of κ in x. torch.linalg.cond computes it from the singular values when p=None, the default, which gives the 2-norm (cond). For optimization, the condition number of the loss's Hessian (λ_max/λ_min) sets how slowly plain gradient descent zig-zags down a long narrow valley; see the linear regression page. Normalizing features and layers mostly works by bringing this number down.

Where they appear in ML

• PCA is the SVD of the centred data matrix. The rows of Vᵀ are the principal directions, and σᵢ² / (n − 1) is the variance along each one. torch.pca_lowrank centres the columns and calls the randomized SVD (pca_lowrank). scikit-learn's PCA with the full solver centres, calls an exact SVD, and sets explained_variance_ = S**2 / (n_samples - 1) (PCA._fit_full). Its randomized solver uses the same Halko algorithm as PyTorch (_randomized_svd).
• LoRA fine-tunes a frozen weight W by learning a low-rank update ΔW = BA with B: [d_out, r] and A: [r, d_in]. It relies on the bet that the change fine-tuning needs has low rank, even though W itself does not.
• Embeddings are low-rank factors. Truncating the SVD of a user × item ratings matrix, or of a word co-occurrence matrix, gives each row and column a short vector, and their dot products rebuild the table. The worked example below does this. Matrix factorization for recommenders and classic word vectors such as LSA start here; see embeddings.
• Initialization. nn.init.orthogonal_ fills a weight with a random orthogonal matrix, so every singular value is 1 and signals keep their size through deep linear stacks (Saxe et al.). It draws a Gaussian matrix, takes its QR factorization, and fixes the signs using R's diagonal so the result is uniformly distributed (orthogonal_, linalg.qr).
• Least squares and the pseudo-inverse. A⁺ = V Σ⁺ Uᵀ inverts each non-zero σᵢ and leaves the zeros at zero, which gives the minimum-norm least-squares solution (pinv).

Worked example

Implementation

import torch

S = torch.tensor([[2., 1.], [1., 2.]])
L, Q = torch.linalg.eigh(S)              # L = tensor([1., 3.]), ascending
Q @ torch.diag(L) @ Q.T                  # rebuilds S

A = torch.tensor([[3., 0.], [4., 5.]])
U, s, Vh = torch.linalg.svd(A)           # s = tensor([6.7082, 2.2361]), descending
U @ torch.diag(s) @ Vh                   # rebuilds A; U and Vh are unique only up to sign
torch.linalg.cond(A)                     # tensor(3.0000)
torch.linalg.matrix_rank(torch.tensor([[1., 2., 3.], [2., 4., 6.], [1., 0., 1.]]))  # tensor(2)

R = torch.tensor([[5., 4, 1, 1], [4, 5, 1, 2], [1, 1, 5, 4], [2, 1, 4, 5], [3, 3, 3, 3]])
U, s, Vh = torch.linalg.svd(R, full_matrices=False)
k = 2
R2 = U[:, :k] @ torch.diag(s[:k]) @ Vh[:k]    # best rank-2 approximation
((R - R2) ** 2).sum(), (s[k:] ** 2).sum()     # equal: Eckart–Young

U, s, V = torch.svd_lowrank(R, q=2)           # randomized; note V, not Vh

W = torch.empty(256, 512)
torch.nn.init.orthogonal_(W)                  # every singular value is 1
torch.linalg.svdvals(W).aminmax()

A one-sided Jacobi SVD from scratch. It applies plane rotations until every pair of columns is orthogonal:

import math

def svd(a, sweeps=60):                   # a: list of rows, len(a) >= len(a[0])
    m, n = len(a), len(a[0])
    cols = [[a[i][j] for i in range(m)] for j in range(n)]
    v = [[float(i == j) for i in range(n)] for j in range(n)]
    dot = lambda x, y: sum(p * q for p, q in zip(x, y))
    for _ in range(sweeps):
        for p in range(n - 1):
            for q in range(p + 1, n):
                alpha, beta, gamma = dot(cols[p], cols[p]), dot(cols[q], cols[q]), dot(cols[p], cols[q])
                if abs(gamma) < 1e-15 * math.sqrt(alpha * beta) or gamma == 0:
                    continue
                zeta = (beta - alpha) / (2 * gamma)
                t = math.copysign(1, zeta) / (abs(zeta) + math.hypot(zeta, 1))
                c = 1 / math.hypot(t, 1); s = c * t
                for vecs in (cols, v):
                    x, y = vecs[p], vecs[q]
                    vecs[p] = [c * i - s * j for i, j in zip(x, y)]
                    vecs[q] = [s * i + c * j for i, j in zip(x, y)]
    sigma = [math.sqrt(dot(c, c)) for c in cols]          # column lengths are the σ's
    u = [[x / s for x in c] for c, s in zip(cols, sigma)]  # assumes every σ > 0
    return u, sigma, v                                   # unsorted; sort by σ to match torch

Tradeoffs

• Don't form AᵀA. Its condition number is κ(A)², so a singular value of 1e-8 becomes an eigenvalue of 1e-16 and disappears into rounding. Computing the SVD of [[1, 2, 3], [2, 4, 6], [1, 0, 1]] through AᵀA reports a third singular value near 5e-8, and so rank 3, where the direct SVD gives rank 2. Use linalg.svd, lstsq or QR instead.
• Cost. A full SVD of an m × n matrix (m ≥ n) costs O(mn²). That is fine for a 4096 × 4096 layer and hopeless for a million × million table. svd_lowrank costs roughly O(mnq) with q ≪ n, but its error depends on how fast the spectrum decays. Use q a little above the rank you want and keep niter at 2 or more.
• Signs and gradients. Singular and eigen vectors are unique only up to sign (or phase), so different hardware can return different ones. Gradients through U, V or eigenvectors blow up when two singular values or eigenvalues are close. Differentiate through svdvals/eigvalsh when you can.
• Truncation throws information away. Eckart–Young minimizes squared error, which counts every entry the same. It is not best for missing entries, which real ratings tables are mostly made of, or for errors that matter more in some places. That is why recommenders fit factors by gradient descent on the observed entries instead.
• Low rank is a bet. LoRA assumes the update has low rank. If the task needs a high-rank change, raising r helps up to a point, and full fine-tuning can still win.

Check yourself

Connections

References

• Eckart and Young. The approximation of one matrix by another of lower rank. Psychometrika 1, 211–218, 1936.
• Trefethen and Bau. Numerical Linear Algebra. SIAM, 1997: lectures 4–5 on the SVD, 12 on conditioning.
• Strang. Introduction to Linear Algebra, 6th ed., chapters 6–7. Wellesley-Cambridge Press, 2023.
• Halko, Martinsson and Tropp. Finding structure with randomness: probabilistic algorithms for constructing approximate matrix decompositions. 2009. arXiv:0909.4061
• Saxe, McClelland and Ganguli. Exact solutions to the nonlinear dynamics of learning in deep linear neural networks. 2013. arXiv:1312.6120
• Hu et al. LoRA: Low-Rank Adaptation of Large Language Models. 2021. arXiv:2106.09685
• PyTorch v2.15.0: linalg.svd, linalg.eigh, svd_lowrank. scikit-learn 1.9.1: PCA.