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.
#include <cstdio>
#include <vector>
#include <cmath>
#include <algorithm>
using Matrix = std::vector<std::vector<double>>;
using Vector = std::vector<double>;
// PA = LU via Gaussian elimination with partial pivoting -- no linear-algebra library.
void luDecompose(Matrix A, Matrix& L, Matrix& U, std::vector<int>& P) {
int n = static_cast<int>(A.size());
U = A;
L.assign(n, Vector(n, 0.0));
P.resize(n);
for (int i = 0; i < n; ++i) P[i] = i;
for (int col = 0; col < n; ++col) {
int pivot = col;
for (int r = col + 1; r < n; ++r)
if (std::fabs(U[r][col]) > std::fabs(U[pivot][col])) pivot = r;
std::swap(U[col], U[pivot]);
std::swap(L[col], L[pivot]);
std::swap(P[col], P[pivot]);
L[col][col] = 1.0;
for (int row = col + 1; row < n; ++row) {
double factor = U[row][col] / U[col][col];
L[row][col] = factor;
for (int c = col; c < n; ++c) U[row][c] -= factor * U[col][c];
}
}
}
Vector forwardSub(const Matrix& L, const Vector& b) {
int n = static_cast<int>(b.size());
Vector y(n, 0.0);
for (int i = 0; i < n; ++i) {
double s = b[i];
for (int j = 0; j < i; ++j) s -= L[i][j] * y[j];
y[i] = s;
}
return y;
}
Vector backSub(const Matrix& U, const Vector& y) {
int n = static_cast<int>(y.size());
Vector x(n, 0.0);
for (int i = n - 1; i >= 0; --i) {
double s = y[i];
for (int j = i + 1; j < n; ++j) s -= U[i][j] * x[j];
x[i] = s / U[i][i];
}
return x;
}
Vector solve(const Matrix& A, const Vector& b) {
Matrix L, U;
std::vector<int> P;
luDecompose(A, L, U, P);
Vector bPermuted(b.size());
for (size_t i = 0; i < b.size(); ++i) bPermuted[i] = b[P[i]];
return backSub(U, forwardSub(L, bPermuted));
}
Matrix inverse2x2(const Matrix& A) {
double a = A[0][0], b = A[0][1], c = A[1][0], d = A[1][1];
double det = a * d - b * c;
return {{d / det, -b / det}, {-c / det, a / det}};
}
Vector matvec(const Matrix& A, const Vector& x) {
Vector out(A.size(), 0.0);
for (size_t i = 0; i < A.size(); ++i)
for (size_t j = 0; j < x.size(); ++j) out[i] += A[i][j] * x[j];
return out;
}
int main() {
Matrix A = {{2., 1.}, {1., 3.}};
Vector b = {5., 10.};
Vector xSlow = matvec(inverse2x2(A), b); // the WRONG way
Vector xFast = solve(A, b); // the RIGHT way
std::printf("x_slow = [%.4f, %.4f]\n", xSlow[0], xSlow[1]);
std::printf("x_fast = [%.4f, %.4f]\n", xFast[0], xFast[1]);
return 0;
}- 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.