Matrix Decompositions: QR, LU, SVD¶
This notebook compares three fundamental matrix decompositions side by side: LU (Gaussian elimination), QR (orthogonal factorisation), and SVD (singular value decomposition). We verify each decomposition by reconstructing the original matrix, explore their properties, and discuss when to use each.
Key functions: np_lu, np_qr, np_svd, np_det, np_diag, np_norm,
np_matmul, np_transpose, np_inv.
load("numerics")$
load("ax-plots")$
Test Matrix¶
We use a well-conditioned $3 \times 3$ symmetric matrix. Symmetric positive-definite matrices arise naturally in least-squares problems, covariance matrices, and finite-element stiffness matrices. The condition number $\kappa(A) = \|A\| \cdot \|A^{-1}\|$ tells us how sensitive linear systems involving $A$ are to perturbations.
/* A well-conditioned symmetric matrix */
A : ndarray(matrix([4, 1, 2], [1, 3, 1], [2, 1, 5]))$
print("Test matrix A:")$
np_to_matrix(A);
/* Condition number: kappa(A) = ||A|| * ||A^{-1}|| */
kappa : np_norm(A) * np_norm(np_inv(A));
print("Condition number kappa(A) =", kappa)$
print("This is well-conditioned (kappa close to 1 is ideal).")$
Test matrix A: matrix([4.0,1.0,2.0],[1.0,3.0,1.0],[2.0,1.0,5.0]) 4.664979478158368 Condition number kappa(A) = 4.664979478158368 This is well-conditioned (kappa close to 1 is ideal).
LU Decomposition¶
The LU decomposition factors $A$ into:
$$PA = LU$$
where $P$ is a permutation matrix, $L$ is lower triangular with ones on the
diagonal, and $U$ is upper triangular. This is the workhorse behind
np_solve -- once we have $L$ and $U$, solving $Ax = b$ reduces to
two triangular solves (forward and back substitution), each $O(n^2)$.
The permutation $P$ ensures numerical stability by performing partial pivoting.
/* LU decomposition: PA = LU */
[L, U, P] : np_lu(A)$
print("L (lower triangular):")$
np_to_matrix(L);
print("U (upper triangular):")$
np_to_matrix(U);
print("P (permutation):")$
np_to_matrix(P);
L (lower triangular): matrix([1.0,0.0,0.0],[0.25,1.0,0.0],[0.5,0.18181818181818182,1.0]) U (upper triangular): matrix([4.0,1.0,2.0],[0.0,2.75,0.5],[0.0,0.0,3.909090909090909]) P (permutation):
/* Verify: PA = LU, so ||PA - LU|| should be ~0 */
lu_error : np_norm(np_sub(np_matmul(P, A), np_matmul(L, U)));
print("Reconstruction error ||PA - LU|| =", lu_error)$
0.0 Reconstruction error ||PA - LU|| = 0.0
QR Decomposition¶
The QR decomposition factors $A$ into:
$$A = QR$$
where $Q$ is an orthogonal matrix ($Q^T Q = I$) and $R$ is upper triangular. QR is the foundation of:
- Least-squares solvers (more numerically stable than the normal equations)
- The QR algorithm for computing eigenvalues
- Gram-Schmidt orthogonalisation
The key property is that $Q$ preserves norms: $\|Qx\| = \|x\|$, which means rounding errors are not amplified.
/* QR decomposition: A = QR */
[Q, R] : np_qr(A)$
print("Q (orthogonal):")$
np_to_matrix(Q);
print("R (upper triangular):")$
np_to_matrix(R);
Q (orthogonal):
matrix([0.8728715609439694,-0.26726124191242434,-0.408248290463863],
[0.2182178902359924,0.9621404708847276,-0.16329931618554513],
[0.4364357804719848,0.053452248382484815,0.8981462390204986])
R (upper triangular):
/* Verify: A = QR */
qr_error : np_norm(np_sub(A, np_matmul(Q, R)));
print("Reconstruction error ||A - QR|| =", qr_error)$
/* Check orthogonality: Q^T Q should be the identity */
QtQ : np_matmul(np_transpose(Q), Q)$
ortho_error : np_norm(np_sub(QtQ, np_eye(3)));
print("Orthogonality error ||Q^T Q - I|| =", ortho_error)$
print("Q^T Q (should be identity):")$
np_to_matrix(QtQ);
6.661338147750939e-16 Reconstruction error ||A - QR|| = 6.661338147750939e-16 2.5051018131878933e-16 Orthogonality error ||Q^T Q - I|| = 2.5051018131878933e-16 Q^T Q (should be identity):
Singular Value Decomposition (SVD)¶
The SVD factors $A$ into:
$$A = U \Sigma V^T$$
where $U$ and $V$ are orthogonal and $\Sigma = \text{diag}(\sigma_1, \sigma_2, \sigma_3)$ contains the singular values in descending order. The SVD is the most informative decomposition:
- $\sigma_1 / \sigma_n$ gives the condition number
- The number of nonzero singular values is the rank
- Truncating small singular values gives the best low-rank approximation
- The columns of $U$ and $V$ are the left and right singular vectors
/* SVD: A = U * diag(S) * Vt */
[U_svd, S, Vt] : np_svd(A)$
print("U (left singular vectors):")$
np_to_matrix(U_svd);
print("Singular values S:", np_to_list(S))$
print("Vt (right singular vectors, transposed):")$
np_to_matrix(Vt);
U (left singular vectors):
matrix([-0.5910090485061037,0.3279852776056822,0.736976229099578],
[-0.3279852776056818,0.7369762290995778,-0.5910090485061041],
[-0.7369762290995784,-0.5910090485061039,-0.32798527760568147])
Singular values S: [7.048917339522306,2.64310413210779,2.3079785283699046]
Vt (right singular vectors, transposed):
/* Reconstruct A from U * diag(S) * Vt */
S_mat : np_diag(np_to_list(S))$
A_reconstructed : np_matmul(np_matmul(U_svd, S_mat), Vt)$
svd_error : np_norm(np_sub(A, A_reconstructed));
print("Reconstruction error ||A - U*diag(S)*Vt|| =", svd_error)$
/* Condition number from singular values */
svals : np_to_list(S)$
kappa_svd : svals[1] / svals[length(svals)];
print("Condition number from SVD (sigma_max / sigma_min) =", kappa_svd)$
2.9394713984828794e-15 Reconstruction error ||A - U*diag(S)*Vt|| = 2.9394713984828794e-15 3.0541520438237635 Condition number from SVD (sigma_max / sigma_min) = 3.0541520438237635
Symbolic Determinant vs Numeric [S+N]¶
The determinant connects to all three decompositions:
- From LU: $\det(A) = \det(P)^{-1} \cdot \prod u_{ii}$
- From SVD: $|\det(A)| = \prod \sigma_i$
We compute the determinant symbolically (exact integer arithmetic) and
numerically (np_det via LAPACK), and verify they agree.
/* Symbolic determinant (exact integer result) */
det_sym : determinant(matrix([4, 1, 2], [1, 3, 1], [2, 1, 5]));
print("Symbolic determinant =", det_sym)$
/* Numeric determinant via LAPACK */
det_num : np_det(A);
print("Numeric determinant =", det_num)$
/* Product of singular values gives |det| */
det_svd : apply("*", np_to_list(S));
print("Product of singular values =", det_svd)$
print("(matches |det| since A is positive definite)")$
43 Symbolic determinant = 43 43.0 Numeric determinant = 43.0 43.00000000000001 Product of singular values = 43.00000000000001 (matches |det| since A is positive definite)
Comparing Decompositions¶
Each decomposition has its strengths:
| Decomposition | Factors | Cost | Best for |
|---|---|---|---|
| LU | $PA = LU$ | $\frac{2}{3}n^3$ | Solving $Ax = b$, computing determinants |
| QR | $A = QR$ | $\frac{4}{3}n^3$ | Least squares, eigenvalue algorithms, stability |
| SVD | $A = U\Sigma V^T$ | $\sim 11n^3$ | Rank, condition number, low-rank approximation |
Rules of thumb:
- For solving linear systems ($Ax = b$): use LU (fastest, what
np_solvedoes) - For least squares ($\min \|Ax - b\|$): use QR (better numerical stability than normal equations)
- For understanding the matrix (rank, condition, principal components): use SVD
- QR and SVD work on non-square matrices; LU requires square matrices
Non-Square Matrix¶
LU decomposition requires a square matrix, but QR and SVD generalise naturally to $m \times n$ matrices. This is essential for least-squares problems where the system is over-determined ($m > n$) or under-determined ($m < n$).
Here we take a $4 \times 2$ matrix (tall, over-determined shape) and show that QR and SVD produce valid factorisations.
/* Non-square matrix: 4x2 (tall) */
B : ndarray(matrix([1, 2], [3, 4], [5, 6], [7, 8]))$
print("B (4x2):")$
np_to_matrix(B);
print("Shape:", np_shape(B))$
B (4x2): matrix([1.0,2.0],[3.0,4.0],[5.0,6.0],[7.0,8.0]) Shape: [4,2]
/* QR of non-square matrix */
[Q_B, R_B] : np_qr(B)$
print("QR decomposition of B:")$
print("Q shape:", np_shape(Q_B))$
np_to_matrix(Q_B);
print("R shape:", np_shape(R_B))$
np_to_matrix(R_B);
qr_B_error : np_norm(np_sub(B, np_matmul(Q_B, R_B)));
print("Reconstruction error ||B - QR|| =", qr_B_error)$
QR decomposition of B:
Q shape: [4,2]
matrix([0.10910894511799629,0.8295150620062542],
[0.3273268353539886,0.4391550328268387],
[0.5455447255899811,0.048795003647425804],
[0.7637626158259735,-0.34156502553198537])
R shape: [2,2]
matrix([9.16515138991168,10.910894511799622],[0.0,0.9759000729485366])
7.483924389696347e-15
Reconstruction error ||B - QR|| = 7.483924389696347e-15
/* SVD of non-square matrix */
[U_B, S_B, Vt_B] : np_svd(B)$
print("SVD of B:")$
print("U shape:", np_shape(U_B))$
np_to_matrix(U_B);
print("Singular values:", np_to_list(S_B))$
print("Vt shape:", np_shape(Vt_B))$
np_to_matrix(Vt_B);
/* Reconstruct: U * diag(S) * Vt */
S_B_mat : np_diag(np_to_list(S_B))$
B_reconstructed : np_matmul(np_matmul(U_B, S_B_mat), Vt_B)$
svd_B_error : np_norm(np_sub(B, B_reconstructed));
print("Reconstruction error ||B - U*diag(S)*Vt|| =", svd_B_error)$
/* Rank from singular values */
print("Matrix rank (number of nonzero singular values) = 2")$
print("The columns of B are linearly independent.")$
SVD of B:
U shape: [4,2]
matrix([-0.15248323331020122,-0.8226474722256603],
[-0.34991837180796403,-0.42137528768457977],
[-0.5473535103057272,-0.020103103143501926],
[-0.7447886488034903,0.3811690813975742])
Singular values: [14.269095499261486,0.6268282324175426]
Vt shape: [2,2]
matrix([-0.6414230279950722,-0.7671873950721771],
[0.7671873950721771,-0.6414230279950722])
8.79251776862914e-15
Reconstruction error ||B - U*diag(S)*Vt|| = 8.79251776862914e-15
Matrix rank (number of nonzero singular values) = 2
The columns of B are linearly independent.
Summary¶
All three decompositions reconstruct their input matrix to machine precision. The choice between them depends on the downstream task:
- LU is the fastest and is used internally by
np_solveandnp_det. It requires a square matrix and uses partial pivoting for stability. - QR is the best choice for least-squares problems and is more numerically stable than LU when the matrix is moderately ill-conditioned.
- SVD is the most expensive but the most informative. It reveals the rank, condition number, and principal directions of the matrix. It works on any shape and is the gold standard for ill-conditioned or rank-deficient problems.