Matrix Decompositions
22 min
Linear Algebra's course taught vectors, matrices, and eigenvalues from
first principles in plain Python — but explicitly stopped short of
matrix DECOMPOSITIONS (LU, QR, SVD), noting they aren't
reasonably hand-implementable the way a dot product or a 2×2
eigenvalue formula is. Real decomposition algorithms need careful
numerical pivoting to stay stable, which is exactly the kind of thing a
tested library like scipy.linalg gets right so you don't have to.
from scipy.linalg import lu
import numpy as np
def lu_reconstructs(matrix):
A = np.array(matrix, dtype=float)
P, L, U = lu(A)
reconstructed = P @ L @ U
return bool(np.allclose(reconstructed, A))
LU decomposition splits a matrix into a permutation P and two triangular factors L and U -- triangular matrices are far cheaper to solve equations with than a general matrix, which is exactly why this decomposition exists: it's the fast path scipy.linalg.solve uses internally for solving linear systems.
Write `lu_reconstructs(matrix)`: use `scipy.linalg.lu` to compute the LU decomposition of `matrix` (`P, L, U = lu(A)`), then return whether `P @ L @ U` reconstructs the original matrix (use `numpy.allclose` to compare, since floating-point decomposition won't be bit-exact).
Linear Algebra's own course deliberately left matrix-decompositions unbuilt, noting LU/QR/SVD 'aren't reasonably hand-implementable the way dot-product/eigenvalues are.' Why not?
You can compute and verify an LU decomposition with scipy.linalg, closing the exact gap Linear Algebra's from-scratch course deliberately left open.