Counting the work: where the n^3 comes from
Across this rung you learned to solve A x = b by first writing A = L U with Gaussian elimination, then sweeping through forward substitution and back substitution. Now we ask the question that decides whether a method is usable on a real problem: how much arithmetic does it actually take? The currency we count in is the flop, a single floating-point operation — one multiply or one add. Counting flops tells you, to leading order, how the running time grows as the matrix gets bigger.
Picture the elimination as a triple-nested loop. To clear the first column you touch nearly all n^2 entries of the matrix; to clear the second column, an (n-1)-by-(n-1) block; the third, smaller still. Summing the work of these shrinking squares — n^2 plus (n-1)^2 plus (n-2)^2 and so on — gives a total of about n^3 / 3 multiply-add pairs, which we report as the order O(n^3). The factorization is the expensive part, full stop. That single number governs everything else you will plan around.
Now contrast that with the two substitution sweeps. Solving a triangular system touches each entry of the triangle once, so forward and back substitution cost only about n^2 flops apiece — a factor of n cheaper than the factorization. For n = 1000 that is the difference between roughly a third of a billion operations to factor and a mere million to substitute. Hold this gap in your mind, because the next idea is built entirely on it.
Factor once, solve many
Here is the payoff of separating A = L U from the substitution. The factors L and U depend only on the matrix A, never on the right-hand side b. So when you must solve A x = b for the same A but a whole stack of different b vectors — which happens constantly, in time-stepping, in optimization, in solving for many load cases of the same structure — you pay the expensive O(n^3) factorization exactly ONCE and then reuse L and U for each new b at the bargain price of O(n^2). This is the principle called factor once, solve many.
factor ONCE: P A = L U cost ~ n^3 / 3 flops
for each b:
solve L y = P b (forward sub) cost ~ n^2 flops
solve U x = y (back sub) cost ~ n^2 flops
ten right-hand sides, n = 1000:
refactor every time : 10 * n^3/3 ~ 3.3 billion flops
factor once : n^3/3 + 10 * 2 n^2 ~ 0.35 billion flopsRecall from the pivoting guide that the factorization you actually store is P A = L U, where P is the partial-pivoting permutation that kept things stable. That changes nothing here: for each new b you simply apply the same permutation first, then run the two cheap substitutions. The recorded pivot order is part of the factorization and gets reused right alongside L and U.
When the matrix is mostly zeros
The O(n^3) cost assumed a DENSE matrix, every entry possibly nonzero. But the giant linear systems that arise in practice — from discretizing a differential equation on a grid, from a circuit, from a road network — are almost always sparse: each row has only a handful of nonzeros, because each unknown is coupled to just a few neighbours. A million-by-million matrix from a 3D grid might have only about five nonzeros per row. Treating it as dense would demand storing 10^12 numbers and doing 10^18 flops — utterly hopeless. Exploiting the zeros is not an optimization here; it is the only way the problem fits in a computer at all.
The cleanest special case is a banded matrix, whose nonzeros all sit on a few diagonals near the main diagonal. The most extreme is tridiagonal — only the main diagonal and its two neighbours — which arises whenever each grid point talks only to its immediate left and right. For tridiagonal systems the Thomas algorithm is just Gaussian elimination with all the guaranteed-zero work skipped, and it solves the whole system in O(n) flops: cost LINEAR in the size, not cubic. A problem that would be hopeless dense becomes trivial.
For a general sparse pattern that is not neatly banded, you reach for a sparse direct solver. It performs the same LU factorization but stores and operates on only the nonzero entries, and it is the routine you should actually call — never write your own. But sparsity hides a subtle trap, and the next section is about that trap.
Fill-in: the price of factoring a sparse matrix
Here is the catch that surprises everyone. Even when A is sparse, its factors L and U usually are NOT. Elimination creates brand-new nonzeros in positions that were zero in A — entries that get filled in as the algorithm combines rows. This is fill-in, and it can be ruinous: a sparse matrix whose factors are nearly dense gives back the O(n^3) cost and the O(n^2) storage you were trying to escape. The sparsity of A is no guarantee of the sparsity of L and U.
The remedy is wonderfully clever: REORDER the unknowns before factoring. Permuting the rows and columns of A leaves the solution unchanged (it just relabels the equations and variables), but it can dramatically change how much fill-in elimination produces. A famous toy example is the 'arrow' matrix, dense along its first row and column and zero elsewhere; eliminate in the given order and it fills in completely, but flip the order so the dense row and column come LAST and there is essentially no fill-in at all. Sparse direct solvers run a fill-reducing ordering heuristic — names like approximate minimum degree and nested dissection — as a cheap planning step before any arithmetic. Choosing a good permutation can be the difference between a solve finishing in a second and never finishing.
Never invert the matrix
Now the title's promised heresy. In a linear-algebra course you write the solution of A x = b as x = A^(-1) b, and it is tempting to compute it that way: form the inverse, multiply. Resist. For solving a linear system, explicitly forming the matrix inverse is slower, less accurate, and wasteful of memory — wrong on all three counts at once. The clean notation x = A^(-1) b is a mathematical statement of WHAT x is, not a recipe for HOW to get it.
Take cost first. Computing A^(-1) means solving A X = I for the n columns of the identity — n right-hand sides — and then you still have to multiply A^(-1) by b. Factoring once and running two substitutions on the single b you care about is several times cheaper and skips the final matrix-vector product entirely. The inverse does strictly more work to reach the same x. And the sparsity story makes it far worse: the inverse of a sparse matrix is almost always DENSE, so inverting a large sparse system manufactures a matrix that will not even fit in memory, while its factors stay sparse.
Accuracy seals the verdict. Forming the inverse and then multiplying introduces extra rounding and is generally less backward stable than solving directly with the factorization — you pay more digits of error for the privilege of doing more work. And neither route can beat the problem's own conditioning: from the conditioning rung, if kappa(A) is around 10^8 you lose about 8 of your ~16 double-precision digits no matter which method you use. The inverse cannot rescue an ill-conditioned problem; it only adds avoidable error on top. So the rule is blunt: to solve A x = b, factor and substitute. The only honest reason to form A^(-1) is the rare case where you genuinely need the inverse's entries themselves — and even then, think twice.
Pulling the rung together
Step back and see what this rung built. You can factor a matrix into triangular pieces, solve triangular systems by substitution, pivot for stability, halve the work with Cholesky when the matrix is symmetric positive-definite, and now reason about cost, exploit sparsity, and refuse the inverse. The thread tying it all together is a single design pattern: do the expensive factorization once, keep its cheap structure, and reuse it. That pattern — and the habit of asking 'what does this cost as n grows?' — will carry straight into the next rung.