LU Decomposition and Singular Value Decomposition

LU decomposition

Gaussian elimination without row swaps factors a square matrix as

$$A = LU,$$

where $U$ is the upper-triangular result of elimination and $L$ is unit lower triangular holding the multipliers: $\ell_{ij}$ is the multiple of row $j$ subtracted from row $i$.

Solving with LU

$A\mathbf x = \mathbf b \iff L(U\mathbf x) = \mathbf b$:

1. Forward substitution: solve $L\mathbf y = \mathbf b$.
2. Back substitution: solve $U\mathbf x = \mathbf y$.

Factoring costs $\frac23n^3$ once; each new right-hand side then costs only $2n^2$. This is why libraries factor once and solve many times.

Pivoting and variants

• If a zero pivot appears, swap rows: $PA = LU$ with a permutation matrix $P$. Every invertible matrix has such a factorisation; partial pivoting also improves numerical stability.
• An LU factorisation without pivoting exists iff all leading principal minors are non-zero.
• LDU: factor the pivots out of $U$: $A = LDU$ with both triangular factors unit-diagonal; for symmetric $A$, $A = LDL^\top$.
• Cholesky: a symmetric positive definite $A$ factors as $A = LL^\top$ (about half the work of LU, no pivoting needed). Its existence is itself a test for positive definiteness.

Singular value decomposition

Every real $m\times n$ matrix factors as

$$A = U\Sigma V^\top = \sum_{i=1}^r \sigma_i\,\mathbf u_i\mathbf v_i^\top,$$

• $U$ ($m\times m$) orthogonal — left singular vectors, eigenvectors of $AA^\top$;
• $V$ ($n\times n$) orthogonal — right singular vectors, eigenvectors of $A^\top A$;
• $\Sigma$ ($m\times n$) diagonal with $\sigma_1 \ge \sigma_2 \ge \dots \ge \sigma_r > 0$, the singular values, $\sigma_i = \sqrt{\lambda_i(A^\top A)}$, and $r = \operatorname{rank}A$.

Computing it by hand

1. Form $A^\top A$; find its eigenvalues $\lambda_i$ and orthonormal eigenvectors $\mathbf v_i$.
2. $\sigma_i = \sqrt{\lambda_i}$.
3. $\mathbf u_i = A\mathbf v_i/\sigma_i$ for each non-zero $\sigma_i$ (complete to an orthonormal basis if needed).

Geometry

Any linear map is a rotation/reflection ($V^\top$), then an axis-aligned stretch ($\Sigma$), then another rotation/reflection ($U$). The unit sphere maps to an ellipsoid whose semi-axes have lengths $\sigma_i$ along $\mathbf u_i$.

What the SVD gives you

Low-rank approximation (Eckart–Young)

The best rank-$k$ approximation to $A$ — in both the spectral and Frobenius norms — is the truncated SVD

$$A_k = \sum_{i=1}^k\sigma_i\mathbf u_i\mathbf v_i^\top, \qquad \lVert A - A_k\rVert_2 = \sigma_{k+1}, \qquad \lVert A - A_k\rVert_F = \sqrt{\sum_{i>k}\sigma_i^2}.$$

Used for image compression, latent semantic analysis, recommender systems (matrix factorisation) and denoising. Storing $A_k$ takes $k(m + n + 1)$ numbers instead of $mn$.

Pseudo-inverse

$$A^+ = V\Sigma^+U^\top, \qquad \Sigma^+ = \operatorname{diag}(1/\sigma_1, \dots, 1/\sigma_r, 0, \dots).$$

$A^+\mathbf b$ is the least-squares solution of $A\mathbf x \approx \mathbf b$ with the smallest norm. If $A$ has full column rank, $A^+ = (A^\top A)^{-1}A^\top$; if $A$ is invertible, $A^+ = A^{-1}$.

SVD and PCA

For a centred data matrix $X$ ($n$ samples × $p$ features), the covariance is $\frac{1}{n-1}X^\top X = V\frac{\Sigma^2}{n-1}V^\top$. So the principal directions are the right singular vectors $\mathbf v_i$, the variance explained by each is $\sigma_i^2/(n-1)$, and the principal component scores are $XV = U\Sigma$. Libraries compute PCA via the SVD of $X$ for numerical stability. See PCA.

Quick revision