The factorization big picture: LU, QR, $Q\Lambda Q^\top$, $X\Lambda X^{-1}$, $U\Sigma V^\top$

Before this lesson: lesson 05 (eigenvalues, diagonalization and Jordan form), lesson 07 (positive definite matrices, Cholesky and the Schur complement) and lesson 08 (the SVD). Time: about 75 minutes, plus the notebook.

---

Why this matters #

Almost every numerical method follows one pattern: factor once, then do cheap operations with the factors. A solver never computes $A^{-1}$; it computes LU or Cholesky and then two triangular solves. Least squares uses QR. PCA uses the SVD or $Q\Lambda Q^\top$. Linear dynamics use $X\Lambda X^{-1}$. Gaussian processes, Kalman filters and quasi-Newton methods (L-BFGS) update inverses with the Woodbury formula. Knowing which factorization fits which problem, and what it costs, is a standard interview question and a daily decision in scientific computing (implicit time stepping, Gaussian processes, data assimilation).

1. Predict first #

2. The picture #

Each factorization writes $A$ as a product of matrices of three easy kinds:

• triangular: $Tx = c$ is solved by substitution in about $n^2$ flops;
• orthogonal: the inverse is the transpose, and lengths are preserved (no error growth);
• diagonal: $n$ independent scalar operations.

Read each factorization as a story about the map $x \mapsto Ax$:

Factorization	Story
$A = LU$	elimination: $U$ is what is left after the row operations, $L$ records them
$A = QR$	an orthonormal basis $Q$ for the columns; $R$ holds the coordinates of the columns in it
$S = Q\Lambda Q^\top$	rotate to the eigenvector axes, scale each axis, rotate back
$A = X\Lambda X^{-1}$	change to a (possibly skewed) eigenvector basis, scale, change back
$A = U\Sigma V^\top$	rotate, scale, rotate, with two different orthonormal bases

The last three rows are easy to mix up. The figure puts them side by side for a $2\times 2$ matrix.

Interactive figure: loading…

3. Definitions and main results #

LU is elimination, written as matrices #

Simplest case, $2\times 2$ with $a_{11} \neq 0$. Subtract $\ell = a_{21}/a_{11}$ times row 1 from row 2:

So $A = E^{-1}U = LU$ with $L = \begin{pmatrix}1 & 0\\ \ell & 1\end{pmatrix}$: inverting $E$ only flips the sign of the multiplier. For $n\times n$, $L$ is unit lower triangular (ones on the diagonal) and holds every multiplier $\ell_{ij}$ below its diagonal; $U$ is upper triangular with the pivots on its diagonal. Cost: $\tfrac23 n^3$ flops for the factorization, then $2n^2$ per right-hand side (forward substitution $Ly = b$, back substitution $Ux = y$).

Interactive figure: loading…

Why the multipliers go into $L$ unchanged #

Take $A = \begin{pmatrix} 1 & 1 & 1 \\ 2 & 3 & 3 \\ 4 & 7 & 8 \end{pmatrix}$ (type it into the figure with Edit A to watch the steps). Elimination uses three multipliers:

1. $\ell_{21} = 2/1 = 2$: row 2 $-$ 2 · row 1 gives $(0, 1, 1)$.
2. $\ell_{31} = 4/1 = 4$: row 3 $-$ 4 · row 1 gives $(0, 3, 4)$.
3. $\ell_{32} = 3/1 = 3$: row 3 $-$ 3 · (new row 2) gives $(0, 0, 1)$.

So $U = \begin{pmatrix} 1 & 1 & 1 \\ 0 & 1 & 1 \\ 0 & 0 & 1 \end{pmatrix}$. Each step is a matrix $E_{ij}$: the identity with $-\ell_{ij}$ in position $(i, j)$. Then $E_{32}E_{31}E_{21}A = U$.

The product of the steps mixes the multipliers.

The $(3,1)$ entry is $2$, not $-\ell_{31} = -4$. It is $-\ell_{31} + \ell_{32}\ell_{21} = -4 + 6$: step 3 subtracts a multiple of the new row 2, which already contains $-2$ times row 1.

The inverse does not mix them.

Why, in one sentence. When row $j$ is used as a pivot row, it is already final: it is row $j$ of $U$. So each row of $A$ is its final row plus the multiples of pivot rows that were subtracted from it:

and this is $A = LU$, read row by row (lesson 01, row view). Check row 3: $(0, 0, 1) + 4\,(1, 1, 1) + 3\,(0, 1, 1) = (4, 7, 8)$ ✓.

Theorem (LU without row exchanges). Let $A$ be invertible, and let $A_k$ be its leading $k\times k$ block. $A = LU$ (unit lower $L$, upper $U$) exists if and only if $\det A_k \neq 0$ for $k = 1, \dots, n$. It is then unique, and the pivots are $u_{kk} = \det A_k / \det A_{k-1}$ (with $\det A_0 = 1$).

Proof of the theorem

1. (⇒) Write $A = LU$ in blocks. Because $L$ is lower and $U$ upper triangular, the leading blocks satisfy $A_k = L_k U_k$, so $\det A_k = u_{11}\cdots u_{kk}$ ($\det L_k = 1$). $A$ is invertible, so $U$ is, so every $u_{ii} \neq 0$ and every $\det A_k \neq 0$.
2. (⇐) Run elimination and assume the first $k-1$ pivots were nonzero. The row operations used so far (row $i$ minus a multiple of an earlier row) do not change $\det A_k$. After them, the leading $k\times k$ block is upper triangular with diagonal $u_{11}, \dots, u_{k-1,k-1}$ and the current entry in position $(k,k)$.
3. So that entry equals $\det A_k/\det A_{k-1} \neq 0$: the $k$-th pivot exists. By induction, elimination never meets a zero pivot, and $u_{kk} = \det A_k/\det A_{k-1}$.
4. Uniqueness: if $L_1U_1 = L_2U_2$, then $L_2^{-1}L_1 = U_2U_1^{-1}$. The left side is unit lower triangular, the right side is upper triangular, so both equal $I$. $\blacksquare$

In the zero pivot case of the figure, $A_2 = \begin{pmatrix}1 & 2\\ 2 & 4\end{pmatrix}$ has $\det A_2 = 0$, so the theorem says the second pivot must be 0.

Partial pivoting. A zero pivot stops elimination even when $A$ is invertible ($\begin{pmatrix}0&1\\1&1\end{pmatrix}$), and a tiny pivot destroys accuracy (section 4). So at step $k$, swap in the row with the largest $|a_{ik}|$, $i \ge k$. Then every multiplier has $|\ell_{ij}| \le 1$, and for every square $A$ there is a permutation matrix $P$ with $PA = LU$ (a column that is already zero is skipped). The worst-case growth of the entries is $2^{n-1}$, but in practice it is small, and LU with partial pivoting is the standard dense solver (np.linalg.solve calls LAPACK gesv).

Interactive figure: loading…

Factor once, solve many: never form $A^{-1}$ #

Once $PA = LU$ is known, $Ax = b$ becomes $LUx = Pb$, and two triangular solves finish it: first $Ly = Pb$ (forward substitution), then $Ux = y$ (back substitution). Each costs about $n^2$ flops. The expensive work is all in the factorization, and it does not depend on $b$.

For the predict question ($n = 1000$, 500 right-hand sides):

Plan	Once	Per right-hand side	Total for 500
factor $PA = LU$, then two triangular solves	$\tfrac23 n^3 \approx 0.67\times 10^9$	$2n^2 = 2\times 10^6$	$\approx 1.7\times 10^9$
form $A^{-1}$, then multiply $A^{-1}b$	$\approx 2n^3 = 2\times 10^9$	$2n^2 = 2\times 10^6$	$\approx 3\times 10^9$

The inverse needs LU plus $n$ more solves (one for each column of $I$), so it is never cheaper. It is also less accurate: $x = A^{-1}b$ usually has a larger residual $\|Ax - b\|$ than the triangular solves (lesson 13). Use scipy.linalg.lu_factor once, then lu_solve for each $b$ (or pass all $b$ as the columns of one matrix).

Cholesky: LU for symmetric positive definite matrices #

For SPD $S$, elimination needs no row exchanges and every pivot is positive (lesson 07). Symmetry gives $U = DL^\top$ with $D = \operatorname{diag}(\text{pivots})$, so $S = LDL^\top = (LD^{1/2})(LD^{1/2})^\top = GG^\top$. Cost $\tfrac13 n^3$: half of LU, because symmetry means only half of the matrix has to be computed.

QR: orthogonalize the columns #

Any $A \in \mathbb{R}^{m\times n}$ with $m \ge n$: $A = QR$, with $Q$ ($m\times n$, orthonormal columns) and $R$ ($n\times n$, upper triangular). Gram–Schmidt gives $r_{ij} = q_i^\top a_j$. With full column rank and $\operatorname{diag}(R) > 0$ it is unique. Key link: $A^\top A = R^\top Q^\top Q R = R^\top R$. So, for full column rank and $\operatorname{diag}(R) > 0$, $R^\top$ is the Cholesky factor $L$ of $A^\top A$, obtained without forming $A^\top A$ (which squares the condition number). Use: least squares ($Rx = Q^\top b$), orthonormal bases, the QR eigenvalue algorithm, the range finder of lesson 09. Lesson 12 computes it stably (Householder).

$Q\Lambda Q^\top$: the spectral decomposition #

Needs $S = S^\top$ (real). $Q$ is orthogonal (eigenvectors), $\Lambda$ is real diagonal. $\Lambda$ is unique up to order. $Q$ is unique up to the sign of each column when the eigenvalues are distinct; inside a repeated eigenspace any orthonormal basis works. Use: PCA, quadratic forms, matrix functions $f(S) = Qf(\Lambda)Q^\top$, normal modes, PD tests.

$X\Lambda X^{-1}$: the eigendecomposition #

Needs $n$ independent eigenvectors (a diagonalizable $A$). It fails for defective matrices such as $\begin{pmatrix}0 & 1\\ 0 & 0\end{pmatrix}$ (lesson 05: Jordan form). $X$ is not orthogonal in general, $\Lambda$ may be complex for real $A$, and each column of $X$ can be rescaled (for a repeated eigenvalue, any basis of its eigenspace works). Use: powers $A^k = X\Lambda^kX^{-1}$, linear recurrences, linear ODEs $e^{At} = Xe^{\Lambda t}X^{-1}$, stability of dynamics. Numerical warning: when $X$ is ill-conditioned ($A$ is close to defective), formulas with $X^{-1}$ amplify errors.

$U\Sigma V^\top$: the SVD #

Any $m\times n$ matrix. $U$, $V$ orthogonal, $\Sigma$ diagonal with $\sigma_1 \ge \dots \ge 0$. The $\sigma_i$ are unique; for distinct positive $\sigma_i$ the pair $(u_i, v_i)$ is unique up to a common sign. Use: rank, best low-rank approximation, pseudoinverse, condition number, PCA, matrix norms (lesson 11).

How they relate #

• Symmetric PSD $S$: $Q\Lambda Q^\top$ is an SVD. Symmetric indefinite: flip the sign of $u_i$ for each negative $\lambda_i$; then $\sigma_i = |\lambda_i|$. More generally $\sigma_i = |\lambda_i|$ for normal matrices ($A^\top A = AA^\top$).
• Always: $\sigma_{\min} \le |\lambda| \le \sigma_{\max}$ for every eigenvalue, and $\prod_i \sigma_i = |\det A| = \prod_i |\lambda_i|$.
• $A^\top A = R^\top R$ (QR ↔ Cholesky) and $A^\top A = V\Sigma^2V^\top$ (SVD ↔ spectral).

The map #

Costs are leading flop counts for $n\times n$ (Golub and Van Loan); the iterative ones are approximate.

	Needs	Factors	Unique?	Cost	Main uses
$PA = LU$	square	$P$ permutation, $L$ unit lower, $U$ upper	yes, for invertible $A$ and a fixed $P$	$\tfrac23 n^3$	$Ax = b$, many $b$, $\det A$
$S = LL^\top$	SPD	$L$ lower, diagonal $> 0$	yes	$\tfrac13 n^3$	SPD solves, sampling, $\log\det$, PD test
$A = QR$	any $m\times n$, $m \ge n$	$Q$ orthonormal columns, $R$ upper	yes, if full column rank and $r_{ii} > 0$	$\tfrac43 n^3$	least squares, orthonormal bases
$S = Q\Lambda Q^\top$	$S = S^\top$	$Q$ orthogonal, $\Lambda$ real	$\Lambda$ up to order; $Q$ up to signs / within eigenspaces	$\approx\tfrac43 n^3$ ($\Lambda$), $\approx 9n^3$ (with $Q$)	PCA, quadratic forms, $f(S)$
$A = X\Lambda X^{-1}$	$n$ independent eigenvectors	$X$ eigenvectors, $\Lambda$ eigenvalues (maybe complex)	$\Lambda$ up to order; $X$ up to column scaling (distinct $\lambda_i$)	$\approx 10n^3$ ($\Lambda$), $\approx 25 n^3$ (Schur form with vectors)	$A^k$, $e^{At}$, dynamics
$A = U\Sigma V^\top$	any $m\times n$	$U$, $V$ orthogonal, $\Sigma \ge 0$	$\Sigma$ yes; vectors up to paired signs (distinct $\sigma_i > 0$)	$\approx\tfrac83 n^3$ ($\Sigma$), $\approx 21n^3$ ($U$, $\Sigma$, $V$)	rank, low rank, $A^+$, $\kappa$

Interactive figure: loading…

Counting parameters #

An $n\times n$ orthogonal matrix has $n^2$ entries, but $Q^\top Q = I$ is a symmetric matrix equation: $n(n+1)/2$ independent scalar constraints. So it has $n^2 - n(n+1)/2 = n(n-1)/2$ free parameters (rotation angles). A triangular factor has $n(n+1)/2$ entries, a unit triangular one $n(n-1)/2$.

Factorization	Count	Total
$LU$	$\frac{n(n-1)}{2} + \frac{n(n+1)}{2}$	$n^2$
$QR$	$\frac{n(n-1)}{2} + \frac{n(n+1)}{2}$	$n^2$
$X\Lambda X^{-1}$	$(n^2 - n) + n$ (each column of $X$ has a free scale)	$n^2$
$U\Sigma V^\top$	$\frac{n(n-1)}{2} + n + \frac{n(n-1)}{2}$	$n^2$
$Q\Lambda Q^\top$, $LL^\top$	$\frac{n(n-1)}{2} + n$, or $\frac{n(n+1)}{2}$	$\frac{n(n+1)}{2}$

Every general factorization accounts for exactly $n^2$ numbers. The symmetric ones give $n(n+1)/2$, which is exactly the number of free entries of a symmetric matrix. A count that does not match is a sign of a wrong assumption.

Low-rank updates: the Sherman–Morrison–Woodbury formula #

Often the matrix changes a little after you have factored it: one new data point, one new observation, one quasi-Newton step. If the change has low rank $k$, you do not need a new $O(n^3)$ factorization.

Theorem (Woodbury). Let $A$ ($n\times n$) and $C$ ($k\times k$) be invertible, $U$ be $n \times k$ and $V$ be $k\times n$. If $C^{-1} + VA^{-1}U$ is invertible, then

Rank 1 (Sherman–Morrison). With $k = 1$, $C = 1$, $U = u$, $V = v^\top$:

Small check. $A = I$, $u = v = (1, 2)$. Then $1 + u^\top u = 6$ and $(I + uu^\top)^{-1} = I - uu^\top/6 = \begin{pmatrix} 5/6 & -1/3 \\ -1/3 & 1/3 \end{pmatrix}$. Directly: $I + uu^\top = \begin{pmatrix} 2 & 2 \\ 2 & 5 \end{pmatrix}$ has determinant $6$ and inverse $\tfrac16\begin{pmatrix} 5 & -2 \\ -2 & 2 \end{pmatrix}$ ✓.

Why it is cheap. If you can apply $A^{-1}$ fast (it is diagonal, or you already have its LU or Cholesky factors), the only new inverse is the $k\times k$ matrix $C^{-1} + VA^{-1}U$ (the capacitance matrix). The cost drops from $O(n^3)$ to about $O(n^2 k + k^3)$, or $O(nk^2 + k^3)$ for a diagonal $A$.

Proof: multiply and simplify

Put $W = C^{-1} + VA^{-1}U$ and multiply $A + UCV$ by the right-hand side:

1. $(A + UCV)(A^{-1} - A^{-1}UW^{-1}VA^{-1}) = I + UCVA^{-1} - UW^{-1}VA^{-1} - UCVA^{-1}UW^{-1}VA^{-1}$.
2. Group the last two terms: $- U\,(I + CVA^{-1}U)\,W^{-1}VA^{-1}$.
3. $I + CVA^{-1}U = C(C^{-1} + VA^{-1}U) = CW$, so the group is $-UCWW^{-1}VA^{-1} = -UCVA^{-1}$.
4. What is left is $I + UCVA^{-1} - UCVA^{-1} = I$. $\blacksquare$

ML use: ridge regression with few samples. Let $X$ be $n \times d$ with $n \ll d$ (for example 100 samples and $10^4$ features). Ridge needs $w = (X^\top X + \lambda I_d)^{-1}X^\top y$, a $d\times d$ solve. Woodbury with $A = \lambda I_d$, $U = X^\top$, $C = I_n$, $V = X$ turns it into an $n \times n$ solve:

1. $(\lambda I_d + X^\top X)^{-1} = \tfrac1\lambda\big(I_d - X^\top(\lambda I_n + XX^\top)^{-1}X\big)$.
2. Multiply by $X^\top y$ and factor out $X^\top(\lambda I_n + XX^\top)^{-1}$ on the left: $w = \tfrac1\lambda X^\top(\lambda I_n + XX^\top)^{-1}\big[(\lambda I_n + XX^\top) - XX^\top\big]y$.
3. The bracket is $\lambda I_n$, so $w = X^\top(XX^\top + \lambda I_n)^{-1}y$.

The cost goes from $O(d^3)$ to $O(n^2 d + n^3)$. The $n\times n$ matrix $XX^\top$ holds the inner products of the samples: this is the kernel form of ridge regression, the starting point of kernel ridge regression and Gaussian processes. The Kalman filter uses Woodbury in the same way: $(P^{-1} + H^\top R^{-1}H)^{-1} = P - PH^\top(HPH^\top + R)^{-1}HP$, which is the posterior covariance of lesson 07 (only an $m \times m$ matrix is inverted, $m$ = number of observations). L-BFGS applies rank-2 updates of this kind to its inverse Hessian estimate.

Woodbury is the Schur complement read two ways

Take the block matrix $M = \begin{pmatrix} A & U \\ V & -C^{-1} \end{pmatrix}$ and compute the top-left block of $M^{-1}$ in two ways (lesson 07: block elimination).

1. Eliminate the bottom-right block first. The Schur complement of $-C^{-1}$ is $A - U(-C^{-1})^{-1}V = A + UCV$, and the top-left block of $M^{-1}$ is its inverse, $(A + UCV)^{-1}$.
2. Eliminate $A$ first. The Schur complement of $A$ is $-C^{-1} - VA^{-1}U = -W$, and the block-inverse formula gives the top-left block $A^{-1} + A^{-1}U(-W)^{-1}VA^{-1}$.
3. Both are the same block of the same matrix, so they are equal. This is the Woodbury formula.

The same identity gives the push-through identity $X^\top(XX^\top + \lambda I)^{-1} = (X^\top X + \lambda I)^{-1}X^\top$, which you can also check directly: multiply both sides by $(X^\top X + \lambda I)$ on the left and by $(XX^\top + \lambda I)$ on the right.

Answers to section 1. (1) Something else: factor $PA = LU$ once ($\tfrac23 n^3$ flops), then each $b$ costs $2n^2$ (two triangular solves). Forming $A^{-1}$ costs about $2n^3$ and gives less accurate solutions. (2) No: defective matrices have no eigendecomposition. Yes: every matrix, of any shape, has an SVD.

4. Worked example #

One symmetric matrix, all factorizations. $S = \begin{pmatrix}2 & 1\\ 1 & 2\end{pmatrix}$.

• LU: $\ell = 1/2$, so $L = \begin{pmatrix}1 & 0\\ 1/2 & 1\end{pmatrix}$, $U = \begin{pmatrix}2 & 1\\ 0 & 3/2\end{pmatrix}$. Pivots $2, 3/2$; $\det S = 3$.
• Cholesky: $D = \operatorname{diag}(2, 3/2)$, so $G = LD^{1/2} = \begin{pmatrix}\sqrt2 & 0\\ 1/\sqrt2 & \sqrt{3/2}\end{pmatrix}$. Check: $GG^\top = \begin{pmatrix}2 & 1\\ 1 & 1/2 + 3/2\end{pmatrix} = S$. ✓
• QR (Gram–Schmidt): $q_1 = (2,1)/\sqrt5$, $r_{11} = \sqrt5$, $r_{12} = q_1^\top(1,2) = 4/\sqrt5$. The rest is $(1,2) - \tfrac45(2,1) = (-\tfrac35, \tfrac65)$, with length $3/\sqrt5 = r_{22}$, so $q_2 = (-1, 2)/\sqrt5$. Check: $R^\top R = \begin{pmatrix}5 & 4\\ 4 & 5\end{pmatrix} = S^\top S$. ✓
• Spectral = SVD: $\lambda = 3, 1$ with $q_1 = (1,1)/\sqrt2$, $q_2 = (1,-1)/\sqrt2$. $S$ is SPD, so $\sigma = 3, 1$ and $U = V = Q$. Also $X\Lambda X^{-1}$ with $X = Q$.
• Sherman–Morrison: $S = I + uu^\top$ with $u = (1, 1)$, so $S^{-1} = I - uu^\top/(1 + 2) = \tfrac13\begin{pmatrix}2 & -1\\ -1 & 2\end{pmatrix}$. Check: $\det S = 3$ and the $2\times 2$ inverse formula give the same matrix. ✓

A non-symmetric matrix. $A = \begin{pmatrix}2 & 1\\ 0 & 1\end{pmatrix}$ (the triangular case of the first figure).

• Eigenvalues $2, 1$, eigenvectors $(1,0)$ and $(1,-1)$: they meet at $45°$, not $90°$. $X = \begin{pmatrix}1 & 1\\ 0 & -1\end{pmatrix} = X^{-1}$, and $X\Lambda X^{-1} = A$. ✓
• Singular values: $A^\top A = \begin{pmatrix}4 & 2\\ 2 & 2\end{pmatrix}$ has eigenvalues $3 \pm \sqrt5$, so $\sigma_{1,2} = \sqrt{3\pm\sqrt5} \approx 2.288, 0.874$, not $|\lambda| = 2, 1$.
• Checks: $\sigma_1\sigma_2 = \sqrt{9 - 5} = 2 = |\det A|$, and $0.874 \le 1, 2 \le 2.288$. ✓

Why pivot. $\begin{pmatrix}10^{-20} & 1\\ 1 & 1\end{pmatrix}x = \begin{pmatrix}1\\2\end{pmatrix}$ has solution $x \approx (1, 1)$. Without a row exchange, step by step in double precision:

1. The multiplier is $\ell = 1/10^{-20} = 10^{20}$.
2. $u_{22} = 1 - 10^{20}$ rounds to $-10^{20}$: the entry $a_{22} = 1$ is lost.
3. The right-hand side becomes $2 - 10^{20}$, which also rounds to $-10^{20}$.
4. Back substitution: $x_2 = (-10^{20})/(-10^{20}) = 1$, then $x_1 = (1 - x_2)/10^{-20} = 0$.

The true $x_1$ is about $1$. With the rows swapped, $\ell = 10^{-20}$ and everything is accurate. This is the right end of the second figure.

5. In code #

The libraries confirm the worked example. The notebook for this lesson has these cells and the exercises.

import numpy as np
import scipy.linalg as sl

S = np.array([[2.0, 1.0], [1.0, 2.0]])

p, l, u = sl.lu(S)                 # SciPy convention: S = p @ l @ u
print(l, u)                        # l = [[1,0],[0.5,1]],  u = [[2,1],[0,1.5]]
G = np.linalg.cholesky(S)          # lower triangular, positive diagonal
print(G)                           # [[1.414 0], [0.707 1.225]]
Q, R = np.linalg.qr(S)             # Householder (LAPACK geqrf): signs may differ
print(np.abs(R))                   # |R| = [[√5, 4/√5], [0, 3/√5]] = [[2.236 1.789] [0 1.342]]
print(np.allclose(R.T @ R, S.T @ S))   # True: Rᵀ is (up to signs) the Cholesky factor of SᵀS
lam, V = np.linalg.eigh(S)         # symmetric: real λ, orthonormal V (LAPACK syevd)
print(lam, np.linalg.svd(S, compute_uv=False))   # [1. 3.]  [3. 1.]  (SPD: σ = λ)

A = np.array([[2.0, 1.0], [0.0, 1.0]])
lam, X = np.linalg.eig(A)          # general (LAPACK geev): X is not orthogonal
print(lam.real, abs(X[:, 0] @ X[:, 1]))          # [2. 1.]  0.707... (= cos 45°)
print(np.linalg.svd(A, compute_uv=False))        # [2.288 0.874]: not |λ|

# The multipliers go into L unchanged (no row exchanges here).
A3 = np.array([[1., 1, 1], [2, 3, 3], [4, 7, 8]])
E21 = np.array([[1., 0, 0], [-2, 1, 0], [0, 0, 1]])
E31 = np.array([[1., 0, 0], [0, 1, 0], [-4, 0, 1]])
E32 = np.array([[1., 0, 0], [0, 1, 0], [0, -3, 1]])
E = E32 @ E31 @ E21
print(E @ A3)                      # U = [[1,1,1],[0,1,1],[0,0,1]]
print(E)                           # [[1,0,0],[-2,1,0],[2,-3,1]]: the (3,1) entry is mixed
print(np.linalg.inv(E))            # L = [[1,0,0],[2,1,0],[4,3,1]]: the multipliers

# Factor once, solve many right-hand sides.
rng = np.random.default_rng(0)
B = rng.standard_normal((300, 300)) + 300 * np.eye(300)
rhs = rng.standard_normal((300, 500))
lu_piv = sl.lu_factor(B)           # PA = LU, 2n³/3 flops, done once
Xs = sl.lu_solve(lu_piv, rhs)      # 2n² flops per right-hand side
print(np.allclose(B @ Xs, rhs))    # True

# Sherman–Morrison and ridge with n << d (Woodbury).
uu = np.array([1.0, 2.0])
print(np.allclose(np.linalg.inv(np.eye(2) + np.outer(uu, uu)), np.eye(2) - np.outer(uu, uu) / 6))  # True
Xd, yd, lam_r = rng.standard_normal((5, 40)), rng.standard_normal(5), 0.3
w_primal = np.linalg.solve(Xd.T @ Xd + lam_r * np.eye(40), Xd.T @ yd)   # 40 x 40 solve
w_dual = Xd.T @ np.linalg.solve(Xd @ Xd.T + lam_r * np.eye(5), yd)       # 5 x 5 solve
print(np.allclose(w_primal, w_dual))   # True

# A tiny pivot without row exchange, done by hand for the 2x2 system.
e = 1e-20
ell = 1 / e                        # multiplier 1e20
u22 = 1 - ell * 1                  # rounds to -1e20: the "1" is lost
x2 = (2 - ell * 1) / u22           # 1.0
x1 = (1 - x2) / e                  # 0.0, but the true x1 is about 1
print(x1, x2, np.linalg.solve([[e, 1], [1, 1]], [1, 2]))   # 0.0 1.0 [1. 1.]

What np.linalg.svd runs

np.linalg.svd calls LAPACK gesdd: Householder bidiagonalization, then a divide-and-conquer solver for the bidiagonal SVD when singular vectors are wanted. Reference LAPACK (bdsdc) solves only the small blocks (size at most 25) with QR iteration, and uses the dqds algorithm when only singular values are wanted. So the main solver for the bidiagonal matrix is divide and conquer, not QR iteration.

6. Where it shows up #

• Implicit time stepping and Newton solvers (PDE models, climate and weather dynamical cores): a sparse LU of the Jacobian is reused for many steps while the Jacobian is frozen. Sparse factorizations also create fill-in, which decides memory.
• Gaussian processes and kriging: Cholesky of $K + \sigma^2 I$ gives both the solve and $\log\det = 2\sum_i \log L_{ii}$. Its $O(n^3)$ cost is the reason for sparse and low-rank approximations; low-rank (Nyström, inducing-point) approximations $K \approx UCU^\top$ are then inverted with the Woodbury formula.
• Uncertainty and data assimilation: correlated perturbations $x = \mu + Lz$ with $z \sim \mathcal{N}(0, I)$ and $\Sigma = LL^\top$; square-root Kalman filters. The Kalman update inverts only an $m\times m$ matrix (Woodbury), with $m$ the number of observations.
• Regression: QR (or SVD when the matrix is rank deficient). With more features than samples, ridge is solved in the kernel form $X^\top(XX^\top + \lambda I)^{-1}y$.
• Optimization: BFGS and L-BFGS update an inverse Hessian estimate by rank-2 corrections, which is the Sherman–Morrison–Woodbury idea.
• Dynamics: eigenvalues of a linearization decide stability; DMD and Koopman methods fit $X\Lambda X^{-1}$ models to data. For non-normal $A$ the transient behavior can differ a lot from what the eigenvalues suggest.
• Deep learning: torch.linalg has all of these and differentiates through them. The gradients of eigh contain factors $1/(\lambda_i - \lambda_j)$ (for svd, $1/(\sigma_i^2 - \sigma_j^2)$), so they become unstable when eigenvalues or singular values are close.

7. Common mistakes #

• Computing $A^{-1}$ to solve $Ax = b$. About 3× the work of LU, and less accurate.
• LU without pivoting on a general matrix.
• Building $L$ from the product of the elimination steps. $E_{32}E_{31}E_{21}$ mixes the multipliers; $L$ is its inverse, which holds them unchanged.
• Calling eig on a symmetric matrix. Use eigh: real eigenvalues, orthonormal eigenvectors, faster.
• Assuming $A = X\Lambda X^{-1}$ exists, or that $X$ is orthogonal.
• Assuming $\sigma_i = |\lambda_i|$ for a non-normal matrix.
• Reading SciPy's LU as $PA = LU$. scipy.linalg.lu returns $A = p\,l\,u$, so $P = p^\top$.
• Using Woodbury when it does not help. It pays off only when $A^{-1}$ is cheap to apply and $k \ll n$. Repeated rank-1 updates also accumulate rounding error; refactor from time to time.

8. Check yourself #

1. Factor $A = \begin{pmatrix}2 & 1 & 1\\ 4 & -6 & 0\\ -2 & 7 & 2\end{pmatrix}$ as $LU$ and give $\det A$.
   Answer

   $L = \begin{pmatrix}1&0&0\\2&1&0\\-1&-1&1\end{pmatrix}$, $U = \begin{pmatrix}2&1&1\\0&-8&-2\\0&0&1\end{pmatrix}$, $\det A = 2\cdot(-8)\cdot 1 = -16$. (This is the *3×3 example* of the LU figure.)
2. Which factorization? (a) 1000 solves with the same SPD covariance; (b) an overdetermined least-squares fit with $\kappa(A) = 10^8$; (c) $A^{100}$ for a diagonalizable $A$ with a well-conditioned $X$; (d) the best rank-5 approximation.
   Answer

   (a) Cholesky once, then two triangular solves per right-hand side. (b) QR (or SVD); not the normal equations, since $\kappa(A^\top A) = 10^{16}$. (c) $X\Lambda^{100}X^{-1}$. (d) Truncated SVD.
3. A $2\times 2$ matrix has eigenvalues $2$ and $1$ and $\sigma_1 = 2.288$. What is $\sigma_2$?
   Answer

   $\sigma_1\sigma_2 = |\det A| = |\lambda_1\lambda_2| = 2$, so $\sigma_2 = 2/2.288 \approx 0.874$.
4. In $PA = LU$ with partial pivoting, why is $|\ell_{ij}| \le 1$, and why does it matter?
   Answer

   The pivot is the largest entry in absolute value in its column (rows $\ge k$), so each multiplier $a_{ik}/a_{kk}$ has absolute value at most 1. Small multipliers keep the entries of $U$ from growing, which keeps rounding errors small (backward stable in practice).
5. Why does $S = Q\Lambda Q^\top$ have $n(n+1)/2$ parameters and not $n^2$?
   Answer

   $Q$ has $n(n-1)/2$ (orthogonality removes $n(n+1)/2$ of its $n^2$ entries) and $\Lambda$ has $n$: total $n(n+1)/2$, the number of free entries of a symmetric matrix.
6. You know $A^{-1}$ for a $1000 \times 1000$ matrix. Then one entry changes: $a_{37}$ grows by $\delta$. How do you get the new inverse in $O(n^2)$?
   Answer

   The change is the rank-1 matrix $\delta\, e_3 e_7^\top$. Sherman–Morrison with $u = \delta e_3$, $v = e_7$ gives $(A + \delta e_3e_7^\top)^{-1} = A^{-1} - \delta\, (A^{-1}e_3)(e_7^\top A^{-1}) / (1 + \delta\, e_7^\top A^{-1} e_3)$: column 3 and row 7 of $A^{-1}$, one outer product, $O(n^2)$ work. It fails if $1 + \delta (A^{-1})_{73} = 0$: then the new matrix is singular.

Recall cardsAnswer in your head first, then reveal. Add the cards to your review deck to see them again just before you would forget.

1. What do $P$, $L$ and $U$ store in $PA = LU$?

   Show answer

2. Why does each multiplier $\ell_{ij}$ go into $L$ unchanged?

   Show answer

3. Why is partial pivoting used?

   Show answer

4. What are the leading costs of LU and of Cholesky for an $n\times n$ matrix, and of each solve after them?

   Show answer

5. Why should a solver never form $A^{-1}$?

   Show answer

6. Which factorization do you use for a least-squares problem, and why not the normal equations?

   Show answer

7. When does $A = X\Lambda X^{-1}$ exist? Give a matrix where it does not.

   Show answer

8. Which factorization exists for every $m\times n$ matrix?

   Show answer

9. State the Woodbury formula.

   Show answer

10. Count the free parameters of $A = QR$ for an $n\times n$ matrix.

    Show answer

Add these cards to my reviewOpen the review →

9. Interview angle #

• "How do you solve $Ax = b$ for many $b$?" — Factor once (LU, or Cholesky if SPD), then $O(n^2)$ per right-hand side. Never form $A^{-1}$.
• "How do you sample from $\mathcal{N}(\mu, \Sigma)$?" — $x = \mu + Lz$ with $\Sigma = LL^\top$; if $\Sigma$ is only PSD, use $Q\Lambda^{1/2}z$ from eigh.
• "Eigendecomposition vs SVD?" — The SVD always exists, with two orthonormal bases and $\sigma \ge 0$; the eigendecomposition may not exist, uses one (skewed) basis, and describes $A^k$.
• "What do LU, Cholesky, QR and the SVD cost?" — $\tfrac23 n^3$, $\tfrac13 n^3$, $\tfrac43 n^3$, and about $21n^3$ with all singular vectors.
• "You have $10^4$ features and 100 samples. How do you fit ridge?" — Use the kernel form $w = X^\top(XX^\top + \lambda I)^{-1}y$ (Woodbury): a $100\times 100$ solve instead of a $10^4\times 10^4$ one.

10. Practice #

The notebook for this lesson runs in Colab. It has five exercises, each with a check cell and a folded solution:

1. Write LU with partial pivoting, $PA = LU$, and compare it with scipy.linalg.lu.
2. Write forward and back substitution, and solve many right-hand sides with one factorization.
3. Write the Cholesky factorization of an SPD matrix.
4. Use the Cholesky factor to draw correlated samples with a given covariance.
5. Count the free parameters of LU, QR, $Q\Lambda Q^\top$, $X\Lambda X^{-1}$ and the SVD.

The notebook also shows the tiny-pivot failure, the worst-case growth $2^{n-1}$ of partial pivoting, and a demo of the Woodbury formula (ridge with $n \ll d$, timed against the $d \times d$ solve).

11. Watch and read #

4. Factorization into A = LUGilbert Strang · MIT 18.06 Linear Algebra, lecture 4 · 48 min

Gaussian Elimination with Partial Pivoting | Lecture 25 | Numerical Methods for EngineeringJeffrey Chasnov · Numerical Methods for Engineering · 5 min

Lecture 2: Multiplying and Factoring MatricesGilbert Strang · MIT 18.065 (2018) · 48 min

Lecture 18: Counting Parameters in SVD, LU, QR, Saddle PointsGilbert Strang · MIT 18.065 (2018) · 49 min

• Gilbert Strang, MIT 18.06, lecture 4: Factorization into A = LU (course page with notes and exercises).
• Driscoll and Braun, Fundamentals of Numerical Computation, free online: §2.4 LU factorization and §2.6 Row pivoting.
• Golub and Van Loan, Matrix Computations, 4th edition: §3.2 and §3.4 (LU, pivoting), §4.2 (Cholesky), §5.2 (QR), §7.5 and §8.3 (QR algorithm), §8.6 (SVD); §2.1.4 (the Sherman–Morrison–Woodbury formula).
• Nick Higham, "What is the Sherman–Morrison–Woodbury formula?": the formula, a derivation, and where it is used.