LU Decomposition & Solving Linear Systems
How Ax = b is actually solved in production, without ever inverting A.
On this page
Beginner: when you need to solve Ax = b for x, computing a full inverse (section 1.6) is overkill — like buying a whole toolbox to turn one screw. LU decomposition splits A into a Lower-triangular matrix L and an Upper-triangular matrix U, so that A = LU. Triangular systems are cheap to solve directly, without ever forming an inverse.
Intermediate: this is literally the "Gaussian elimination" method taught in introductory algebra, just organized and named so the steps can be reused. Once you have A = LU, solving Ax = b becomes two easy steps: solve Ly = b for y (forward substitution, top to bottom), then solve Ux = y for x (back substitution, bottom to top).
Advanced: plain LU decomposition can fail or become numerically unstable if a "pivot" (a diagonal entry used during elimination) is zero or very small. The fix, partial pivoting, reorders rows during the process, giving PA = LU for some permutation matrix P — this is what every real library actually computes, silently, whenever you call a "solve" function.
L has 1s on its diagonal and zeros above; U has arbitrary values on and above its diagonal and zeros below; P records any row swaps needed for numerical stability.
"Subtract m times row j from row i" is itself a linear operation on the rows of A — which means it can be written as left-multiplying A by an elementary matrix E: the identity matrix (section 1.6) with a single extra −m placed in position (i, j). One elimination step is A → EA.
Elementary matrices are trivially invertible — undoing "subtract m times row j from row i" is just "add m times row j back to row i," so E⁻¹ is E with that one entry negated back to +m.
Eliminate every entry below the diagonal with a sequence E_k⋯E_1, arriving at the upper-triangular result: E_k⋯E_1 A = U. Solving for A:
L is a product of lower-triangular matrices (each E⁻¹ is lower-triangular with 1s on the diagonal), and the product of lower-triangular matrices is always lower-triangular — which is exactly why L comes out lower-triangular, matching the definition above.
Where this is used: this is precisely what a numerical library does internally when you call a solver — Gaussian elimination isn't a separate algorithm from LU decomposition, it is LU decomposition, just narrated step by step instead of packaged into L and U at the end.
Solving 2x + y = 5, x + 3y = 10 by elimination: subtract ½ of row 1 from row 2 to zero out the x-coefficient, giving the "U" row 2.5y = 7.5 → y = 3, then back-substitute into row 1: 2x + 3 = 5 → x = 1. The multiplier you used (½) is exactly the entry that goes into L. This tiny manual example is LU decomposition — the algorithm just formalizes and records every step.
Factor once with lu_factor, then reuse it with lu_solve for as many right-hand sides as you need — this is exactly the performance trick used in iterative simulations that repeatedly solve against the same system matrix.
import numpy as np
from scipy.linalg import lu, lu_factor, lu_solve
A = np.array([[2., 1.], [1., 3.]])
b = np.array([5., 10.])
# The WRONG way (slow, numerically worse at scale):
x_slow = np.linalg.inv(A) @ b
# The RIGHT way — let LAPACK pick the best method (LU under the hood):
x_fast = np.linalg.solve(A, b)
# Explicit LU, if you need to solve against many different b's efficiently:
lu_piv = lu_factor(A)
x_reused = lu_solve(lu_piv, b)
print(x_slow, x_fast, x_reused) # all equal: [1. 3.]- Structural and circuit simulation software solves enormous sparse linear systems via specialized sparse LU decompositions, over and over, as a simulation steps through time.
- Any time a codebase calls
np.linalg.solve,scipy.sparse.linalg.spsolve, or similar, LU decomposition (with pivoting) is almost certainly running underneath.
- Computing
np.linalg.inv(A) @ binstead ofnp.linalg.solve(A, b)— it's slower, less accurate, and considered bad practice in every numerical computing community for exactly the same reason section 1.6 flagged it. - Assuming LU decomposition always exists without pivoting — it can fail on a perfectly invertible matrix if you're unlucky with row order; always use the pivoted version.
Going deeper
LU decomposition costs roughly O(n³/3) operations — about half the cost of computing a full matrix inverse, which is one of the concrete reasons "solve, don't invert" is not just stylistic advice but a real performance and stability rule.
At the master level: for the special case of a symmetric positive-definite matrix (section 1.15), Cholesky decomposition is an even cheaper, more stable variant of LU that exploits the symmetry to do roughly half the work again — this is why covariance-matrix heavy code (Gaussian processes, Kalman filters, Bayesian inference) almost always reaches for Cholesky specifically rather than generic LU.