The factor is not sparse

A sparse matrix is not a compressed dense matrix. Its zeros are structural — they are where they are because the problem has no coupling there — and a matrix arising from a physical discretisation has a fixed number of entries per row however large it gets. The five-point Laplacian on a k×k grid has k² variables and about 5k² nonzeros, forever.

That is the good news, and it is the reason problems with millions of unknowns are posed at all. Storing such a matrix costs O(n). Multiplying by it costs O(n).

Solving with it does not, and the reason is the single fact this field exists to state: eliminating a variable couples every variable it touched to every other variable it touched, and each of those couplings is an entry that was zero in A and is not zero in L.

The rule, in one sentence #

Elimination of variable k adds an edge between every pair of remaining variables that were both connected to k.

That is the whole of it. Think of the matrix as a graph: a vertex per variable, an edge wherever there is a nonzero off-diagonal entry. Eliminating a vertex removes it and turns its neighbourhood into a clique. Every edge that clique adds is a fill entry.

The consequence is immediate and it is what makes the subject hard. A vertex with d neighbours creates up to d(d−1)/2 new edges, so eliminating a single well-connected variable can add more entries than the matrix contained. And those new edges raise the degrees of the vertices they join, so the next elimination is more expensive than it would have been. Fill compounds.

Counting it without doing any arithmetic #

The rule above mentions no numbers, and that is not an accident of the exposition — it is a property of the algorithm, and it is what makes sparse direct solving practical.

The fill pattern can be computed from the graph alone: no additions, no multiplications, no rounding, nothing that depends on the values of the entries. A real sparse solver does exactly this first, in a phase called symbolic factorisation, and it is why such a solver can allocate precisely the right amount of memory before it touches a number.

The site computes it both ways and requires them to agree.

Symbolically, by the clique rule, on a boolean matrix. Numerically, by running Cholesky and counting what came out nonzero. On the 12×12 grid both give 1,739 entries in the factor, against 408 in the matrix’s lower triangle. Not close to each other — the same integer, and the assertion is written as an equality of integers rather than as a tolerance, because an integer is the one kind of measurement on this site that cannot be approximately right.

That agreement is a genuine second route in the sense two routes to a number means. The two computations share nothing but the matrix: one is a walk over a graph, the other is a floating-point factorisation.

Where the two routes could legitimately differ #

Worth stating, because a check that can never fail is not a check and this one has a known escape.

The symbolic count is an upper bound. It asks where a nonzero could appear; the numeric count asks where one did. If a fill entry happens to be exactly zero because the contributing terms cancelled, the numeric count is lower and the two disagree without anything being wrong.

On the model problems that does not happen, and the check is written as an equality for that reason. On a matrix with structural symmetries it can, and a solver that allocated from the symbolic count would simply have reserved slightly too much — which is the right direction for an error in a memory estimate to run.

Why the count is about Cholesky and not about LU #

The essay has been quietly assuming a symmetric positive definite matrix, and the assumption is doing more than tidiness.

For a symmetric positive definite matrix, Cholesky requires no pivoting. That is a theorem, and it is what makes the symbolic phase possible: the elimination order can be chosen entirely on structural grounds, before any value is examined, because no numerical consideration will ever demand a different one. Structure and stability are decoupled.

For a general matrix they are not, and the field becomes considerably harder. The swap that is not optional establishes that elimination without a row interchange fails on matrices it has no business failing on, so a sparse LU has to pivot for stability — and the pivot it wants for stability is rarely the one the ordering wanted for sparsity. A row swap chosen at step k changes the fill pattern of every subsequent step, which means the symbolic prediction is no longer a prediction.

The practical resolutions are all compromises. Threshold pivoting accepts any pivot within a factor of the largest, using the freedom to prefer the sparse choice. Static pivoting commits to the structural ordering and perturbs a pivot that turns out to be too small, then repairs the damage with iterative refinement — which is buying the accuracy back doing structural work rather than precision work.

So the clean two-route agreement in this essay is a symmetric-positive-definite result. It is stated here in the case where it is exactly true, and the general case pays for its generality with a symbolic phase that can only bound what the numeric one will do.

What it costs, fitted rather than quoted #

The usual way to state this is a complexity class, and this site’s habit is to measure the exponent.

Across grid sizes from 5×5 to 12×12, fitting the logarithm of the nonzero count against the logarithm of the dimension:

• the matrix grows with an exponent of 1.04 — linear, as promised, since the stencil is fixed;
• its factor grows with an exponent of 1.49.

The second number is the field. $n^{1.5}$ for a two-dimensional problem is the textbook result and here it is a fitted slope over eight points rather than a citation. The assertions bracket it: above 1.05, so the factor genuinely grows faster than the matrix, and below 2, so it is genuinely better than dense.

The practical translation is unpleasant. A grid ten times finer in each direction has a hundred times the variables and about a thousand times the factor. Memory, not time, is what stops a sparse direct solver, and it stops it suddenly.

“No ordering closes it” is the right verdict and it undersells what an ordering does, which is worth separating because the two facts pull in opposite directions for a caller.

The exponent falls from 1.49 to 1.31 and the gap does not close. Across the natural ordering, reverse Cuthill–McKee, minimum degree and nested dissection the factor’s exponent runs 1.49, 1.41, 1.32 and 1.31, while the matrix’s stays at 1.04–1.05 throughout. So an ordering removes about two-thirds of the excess exponent — from 0.44 above the matrix down to 0.26 — and none of it removes the last third, which is the sense in which the fill is structural rather than incidental.

The counts at one size tell a different story from the exponents, and both are needed. At the 12×12 grid the factor holds 1,354 entries under reverse Cuthill–McKee, 1,026 under minimum degree and 1,413 under nested dissection, against the matrix’s 408 — so minimum degree wins at this size while having the worse exponent of the two good orderings. That is the ordinary shape of a crossover: nested dissection is paying a constant to buy a slope, and at 144 unknowns it has not finished paying. A benchmark run at one grid size would rank them exactly backwards for the sizes where the choice matters.

Why this is the reason iterative methods exist #

The comparison that decides how large problems get solved.

A conjugate gradient iteration costs one matrix-vector product, which for a sparse matrix is O(n). The rate the condition number predicts measures the count at forty for the model problem at n = 40. So the iterative route costs O(kn), with k in the tens or hundreds and — importantly — storage of a fixed handful of vectors.

The direct route costs a factorisation whose storage is the number this essay measures. At n = 144 that is 1,739 entries against the matrix’s 408: a factor of four, which is nothing. At n = 10⁶ on a three-dimensional problem it is the difference between a machine that can hold the problem and one that cannot.

So the decision between iterating and factorising is a memory decision before it is a speed decision, and it is made from the fill count, which is available in advance from the symbolic phase. That is a much more definite basis than the iteration count, which the rate the condition number predicts shows can only be bounded to within an order of magnitude.

The direct method’s compensating virtue #

It would be one-sided to leave it there, because sparse direct solvers are not a legacy technology and there are cases where nothing else will do.

A factorisation is reusable. The expensive phase depends only on the matrix, so a system solved against a thousand right-hand sides pays for the factor once and then pays a triangular solve each time. An iterative method pays in full every time.

It has no convergence question. There is no parameter to tune, no preconditioner to design, no stagnation to diagnose. The answer arrives, backward stable, with the residual the site prints beside every factorisation.

And it does not care about conditioning in the way an iteration does. A κ of 10¹⁰ makes an iterative method very slow and makes elimination merely inaccurate in the way the condition number is an amplifier describes, which is a bound nothing can improve on anyway.

The two approaches meet in the middle, at incomplete factorisations: a factor computed with the fill deliberately discarded, used not as a solver but as a preconditioner. That is changing the condition number on purpose, and it is exactly this essay’s quantity being traded against that essay’s.

What a sparse format actually stores #

The figures on this site draw dense arrays, because at the sizes a picture can show, a real sparse format would be machinery with nothing to teach. What is simulated honestly is the count, and it is worth one paragraph on what that count corresponds to.

Compressed sparse column, the usual format, stores three arrays: the nonzero values, the row index of each, and a pointer per column saying where that column’s entries begin. So an entry costs a double and an integer, and a column costs one more integer. The nnz this essay counts is the length of the first two.

Two consequences follow, and both are why sparse solvers are written the way they are. The pattern must be known before the values are stored, since inserting an entry into the middle of a packed array means moving everything after it — which is precisely why the symbolic phase exists rather than being an optimisation. And the indirection costs: every access is a lookup through the index array, so a sparse matrix-vector product runs at a fraction of the throughput a dense one of the same operation count would, because the hardware cannot predict where the next value is.

That second point is a reason iterative methods are less dominant than the operation counts suggest. An iteration’s cost is one sparse product, and a sparse product is memory-bound rather than compute-bound — while the dense triangular solves inside a direct method run at full speed. Counting floating-point operations understates the direct method’s case, sometimes by a large factor.

What is asserted here #

The two fill counts agree exactly, on five different structures, as integers.

The factor is denser than the matrix, at every position of the ordering slider — which is the claim that there is a subject here at all.

And the factorisation is exact, ‖A − LLᵀ‖/‖A‖ below 10⁻¹³ whichever ordering was used. That one is on the badge and it is doing real work: the entire next essay is about orderings that differ by a factor of six in cost, and the only thing that makes that comparison interesting is that they are all equally correct.

The matrix grows linearly and its factor does not, as fitted exponents with brackets on both.

The refusals: the claim that fill is the same under any two orderings must throw, and so must the claim that a good ordering eliminates fill rather than limiting it. Both do.

Three dimensions, where the argument turns #

Everything measured here is two-dimensional, and the exponent changes in a way that decides how the subject is practised.

For a two-dimensional grid the factor grows as $n^{1.5}$ and the work as $n^2$ — bad, and survivable. For a three-dimensional grid with the same nested-dissection reasoning the factor grows as $n^{4/3}$ and the work as $n^2$, which sounds comparable and is not, because n itself is now k³ rather than k². A modest 1,000³ grid is 10⁹ unknowns, and a factor with $n^{4/3}$ entries is 10¹² numbers.

That is the crossing point for the whole field. Two-dimensional problems are routinely solved directly, and the memory is affordable. Three-dimensional problems of any resolution are not, and this is why essentially every large-scale simulation code — fluid dynamics, electromagnetics, structural analysis at scale — is iterative. Not because iteration is more accurate, and not because it is faster per unknown, but because the factor does not fit.

The exponent in the figure is fitted on two-dimensional grids and says so. Extending it to three was long left undone on the grounds that a spy plot of a thousand-variable matrix is a grey square — which is a good reason not to draw one and no reason at all not to count. A fill count needs no picture, and counting changes the paragraph above.

Minimum degree on both, with symbolicCholesky on the permuted matrix:

Fitted, the exponents are 1.306 in two dimensions and 1.491 in three.

At these sizes the three-dimensional exponent is the larger one, which is the opposite ordering to the asymptotic 1.5 against 4/3 quoted above. Both measured exponents also sit on the wrong side of their own asymptote — below it in 2D, above it in 3D. Neither is a contradiction, and the reason is worth stating: the quoted figures are nested dissection as n → ∞, these are minimum degree at n ≤ 512, and the separator structure the asymptotics rest on is barely present in an 8×8×8 grid. An asymptotic exponent is a statement about a limit and this is a statement about the sizes that fit in a measurement.

What it does mean is that the careful argument above — that n^(4/3) “sounds comparable and is not, because n itself is now k³” — is defending a position that does not need defending at the sizes a practitioner meets. At equal numbers of unknowns the three-dimensional factor is already 1.26 to 1.55 times as full, and the ratio is climbing: 2.04 against 2.57 at n ≈ 64, 2.74 against 4.26 at n ≈ 200.

So three dimensions is worse twice over, and the two reasons are independent. The factor is fuller per unknown — that is this essay’s own rule about coupling, applied to a graph where every vertex has six neighbours instead of four — and there are more unknowns for a given resolution. Only the second was in the paragraph above, and it is the weaker of the two.

assertTheThreeDimensionalFillIsWorseAtEqualN builds both families, fits both exponents, and requires the three-dimensional factor to be the fuller one at every comparable size.

The second one is the one people believe #

That second refusal is worth drawing out, because it catches a misreading that an ordering that does not wait for the numbers could otherwise create.

Shown the arrowhead matrix — where one ordering gives a completely dense factor and another gives no fill at all — the natural conclusion is that fill is an artefact of bad ordering, and a good ordering removes it. It does not. On the 12×12 grid, minimum degree is the best of the four orderings drawn and its factor still holds 1,026 entries against the matrix’s 408. Nested dissection, which has the better asymptotics, gives 1,413.

Every ordering fills in. The best available ordering fills in. The question the field answers is not how to avoid fill but how much of it is unavoidable, and finding the ordering that minimises it is NP-hard, so even that question is answered heuristically.

Which makes the arrowhead a teaching example rather than a representative one, and the essay on ordering opens with it for exactly that reason — and then has to spend most of its length walking the conclusion back.

What the count is not #

One clarification, because “the factor is denser than the matrix” invites a conclusion about accuracy that is not there.

Fill costs memory and it costs arithmetic. It does not cost accuracy. A factorisation under the worst ordering on this page and one under the best are both backward stable, both reconstruct the matrix to about 10⁻¹⁶, and both return the same solution to within rounding. The badge on the hero figure carries that number at every position of the slider precisely so the comparison in the next essay cannot be misread as a comparison of quality.

That separation — a decision that costs a great deal and affects no answer — is unusual on this site. Most of what is measured here is a trade between cost and accuracy: precision against speed, incomplete factorisation against exact, sketching against a full pass. Fill is a pure cost with no accuracy dimension at all, which is what makes ordering the one decision in the subject that can be made on structural grounds alone and be certain of losing nothing.

A function of a sparse matrix, which is worse #

A factor of a sparse matrix fills in. A function of one fills in immediately and completely: e^A has a nonzero wherever there is a path of any length, which on a connected graph is everywhere.