What the symbolic phase can only bound

The single most useful property a sparse direct solver has is that it knows how much memory it needs before it allocates any. The factor is not sparse is where that comes from: eliminating a variable turns its neighbourhood into a clique, every edge that clique adds is an entry of L that was a zero of A, and the whole count follows from the graph with no arithmetic at all — no additions, no multiplications, no rounding.

The site checks it as an integer identity on five structures. Symbolic and numeric agree exactly, and on the matrix used in this field with an ordinary diagonal the numbers are 233 and 233.

That property does not survive pivoting, and this essay is about what does.

Why exactness goes #

The symbolic algorithm needs one thing: the order in which variables will be eliminated. Given that, it can walk the elimination graph and produce the exact pattern.

Under partial or threshold pivoting the order is decided by the values. Row 7 is swapped to position 2 because its entry in column 2 happened to be larger than row 2’s, and whether that happens is a question about the numbers, which the analysis phase does not look at. On the 5×5 grid the pivoted factorisation produces 242 entries where the unpivoted symbolic count says 233, and the assertion that the first predicts the second is fed to the build and refused.

Nine entries out of 233 sounds like a rounding of the problem rather than a change of kind. It is a change of kind, because a solver cannot allocate 233 and hope: an allocation that is nine entries too small fails at the moment the 234th entry appears, which is halfway through the numeric phase, on the largest problem the machine will run that day.

So the analysis phase has to produce something that is guaranteed, and a prediction that is usually right is not it.

The bound, and where it comes from #

The guarantee exists, it is provable, and it is computable from the pattern alone.

George and Ng’s result: for any row permutation P, the sparsity pattern of the U factor of PA is contained in the pattern of the Cholesky factor of AᵀA.

The argument is worth a paragraph because it explains why AᵀA and not something else. Consider the QR factorisation of A. Its R factor is, up to signs, the Cholesky factor of AᵀA — the same object the normal equations produce, which is the road that squares the problem’s subject in a different context. And R’s pattern is invariant to row permutations, because permuting the rows of A does not change AᵀA at all. So R is a structure that every row ordering shares. The theorem is then that LU’s U fits inside it.

Computationally that is three steps and none of them touches a value: form the pattern of AᵀA by asking, for each pair of columns, whether any row has a nonzero in both; run the symbolic Cholesky on that pattern; count. The site forms the pattern from the pattern rather than by multiplying AᵀA in floating point, deliberately — forming the product numerically would let a cancellation decide a structural question, and a structural question decided by rounding is exactly the kind of thing this site refuses to let happen quietly.

Checking it as a bound #

A bound that is never checked is decoration. This one is checked at four grid sizes and at eight strategies each — the seven thresholds plus no pivoting at all — and it holds every time.

Two readings of that table, and the second is the one worth carrying away.

The bound holds. Thirty-two separate factorisations, none of them exceeding it. That is the theorem doing what a theorem does, and it is checked rather than cited.

The bound gets looser as the problem grows. 1.39 at sixteen unknowns and 1.68 at forty-nine, and the direction is monotone across every size drawn. A bound that tightens with size would be an allocation strategy; a bound that loosens is a guarantee and an over-allocation, and the ratio is what a code is paying for the guarantee.

The reason is straightforward once stated: the bound covers every row permutation, and the number of row permutations grows with the matrix while the number that any actual pivoting strategy will produce does not. So the union it takes is over a rapidly growing set of possibilities, of which the elimination realises one.

What codes do about it #

Three strategies, and the third is the one that removed the problem for most users.

Allocate the bound. Guaranteed, simple, and pays the 1.4 to 1.7 times over-allocation. Reasonable when memory is not the binding constraint, which is exactly the situation sparse direct solvers are usually not in.

Allocate dynamically. Grow the data structures as fill appears. This is what most general-purpose unsymmetric sparse solvers do, and it is why their inner loops are more complicated than the symmetric ones: an elimination that may reallocate at any step cannot use the tight indexed arithmetic a fixed structure permits. The cost is not the reallocation, it is the indirection everywhere else.

Fix the pivots in advance. Static pivoting, as structure and stability stop being separable describes: choose the sequence up front from the pattern and the magnitudes as they stand, factorise into exactly that structure whatever the numbers do, and replace any pivot that comes out unusably small with something like ε‖A‖. The analysis phase is then exact again — because the order really is decided in advance — and the price is a factorisation of a matrix that is not quite the one that was asked about.

The third looks like cheating and is the one this site finds most interesting, because it is the same move the mixed-precision field makes. The perturbation introduced by replacing a pivot is a backward error of known size, and iterative refinement against accurately computed residuals converges to the answer of the original problem provided κ times that perturbation stays under one. Which is exactly the condition buying the accuracy back measures, with the cheap factorisation coming from a fixed pivot order instead of from a short mantissa.

The two-currency version, again #

The pattern this whole field keeps producing is worth stating in its general form, since three essays have now arrived at it from different directions.

An exact symbolic phase requires a fixed elimination order. A fixed elimination order requires accepting whatever pivots it produces. Accepting whatever pivots it produces means a backward error that is not bounded in advance. So the exactness of the analysis and the boundedness of the backward error cannot both be had, and every design in this field is a position on that trade:

Dynamic pivoting keeps the backward error small and gives up the exact analysis. Static pivoting keeps the exact analysis and gives up the a priori backward error bound, then buys it back afterwards with refinement. The George–Ng bound is the middle position — keep both, and pay in memory for the guarantee.

None of the three is wrong. What would be wrong is a code that claimed both without paying for either, and the way to tell is the same in every case: look at the residual it returns.

Why the bound is the shape it is, and what would make it tighter #

The bound is a union over row permutations, and a natural question is whether a smaller union would do. Two answers, and the second is where the field’s remaining work is.

A smaller union over the same set is not available. The theorem is tight in the sense that matters: for each entry the bound predicts, there exists some row permutation that fills it. So no argument about this set of permutations gives a smaller structure — the looseness is not slack in the proof, it is the difference between “some permutation fills this” and “the permutation this elimination will choose fills this”.

A smaller set is available, and that is where practical work goes. Once a threshold τ is fixed, not every row permutation is reachable: a row whose candidate is below τ times the column maximum can never be chosen, and which rows those are depends on the values. Codes that want a tighter bound therefore do a partial analysis — one pass over the values to establish which candidates are plausible, then a symbolic phase over the reduced set. That is no longer a bound computed from the pattern alone, which is the property this essay has been about, and it is a good deal tighter.

The trade is the familiar one in a new place. The pattern-only bound is reusable across every matrix with this structure and is loose by 1.7 times; the value-aware bound is tight and has to be recomputed whenever the values change. A time-stepping code with a fixed structure and changing coefficients wants the first; a code solving one matrix once wants the second.

There is a third position, and it is what several modern codes actually do: compute the loose bound, start the factorisation, and reduce the allocation once the elimination has committed to enough pivots that the remaining possibilities are few. That is neither a symbolic phase nor a numeric one and it is why the clean three-phase description of a sparse solve is a description of the symmetric case.

The number nobody prints #

One closing observation, and it is the sort of gap this site exists to notice.

Every sparse solver reports its fill. Almost none of them report the bound they could have computed, so the over-allocation ratio — the 1.39, 1.53, 1.62, 1.68 above — is a quantity that exists inside every analysis phase in the world and is essentially never surfaced.

It is a useful number. It says how much of the allocation was insurance, which is directly what a user planning a larger run wants to know, and it says how much the pivoting strategy is departing from the structure the pattern alone predicts. A run with an over-allocation ratio near one has pivoted almost not at all, which is worth knowing about a matrix. A run near two has pivoted a great deal, which is a hint that the matrix is far from diagonally dominant and that the residual is worth looking at.

The site prints it on the figure because there was no reason not to, and the trajectory — monotone upward across four sizes — is a finding rather than a decoration. It is the kind of measurement that turns a guarantee into a planning basis, and the reason it is missing everywhere else appears to be that nobody counted.

What is asserted here #

The bound holds at every threshold and every grid size, including the unpivoted case — thirty-two checks, each one an inequality on an integer.

And it is loose, by more than 1.2 times at every size, asserted in that direction because a bound reported without its slack invites being read as a prediction.

The unpivoted symbolic count is exact, checked as an integer identity on the same grid with an ordinary diagonal so that the arithmetic is benign and any disagreement would be structural rather than a cancellation. 233 and 233.

The refusals: that the symbolic phase predicts the pivoted fill, which is fed the 233-against-242 case and refuses; and that a bound computed for one matrix holds for another, which is fed the 5×5 grid’s bound applied to the 6×6 grid and refuses. The second is there because a bound that is never tested against the wrong matrix is a number a reader can mistake for a universal constant, and this one is a property of a structure.

What the analysis phase can and cannot say, in one paragraph each #

It can say the factor will fit. That is the guarantee, it is provable, it is computable from the pattern in time proportional to the elimination graph rather than to the arithmetic, and it is what lets a solver refuse a problem before spending an hour on it. A code that reports “this will need 605 entries” and then uses 360 has told the truth.

It cannot say what the factorisation will cost. The operation count depends on the fill that actually occurs, which depends on the pivots, which depend on the values. The bound over-estimates the entries by 1.7 times and therefore over-estimates the arithmetic by considerably more, since the work at each step goes as the square of the degree. So a planning number taken from the symbolic phase is safe for memory and misleading for time.

That split is the practical content of this essay and it is the kind of thing that is obvious once stated and invisible otherwise. The two things a user wants to know before running a large sparse solve are whether it will fit and how long it will take, and the analysis phase answers the first exactly and the second not at all.

A last comparison, to the bound this site already had #

The elimination field carries a bound of the same shape and it is worth putting the two side by side, because the pairing says something about what bounds in this subject are for.

The growth factor’s bound is 2ⁿ⁻¹, which at n = 40 is 5.5·10¹¹, and the bound that is never attained measures a median of 3.23 against it. That is loose by eleven orders of magnitude. The structural bound here is loose by 1.7 times.

Both are correct, both are checked, and they are useful in completely different ways. A bound loose by eleven decades cannot be allocated from, cannot be planned with and cannot be quoted to a user; what it is good for is the existence statement — elimination with partial pivoting is backward stable, with a constant that does not depend on the matrix — and for identifying the matrices that attain it, which turn out to be constructions rather than anything that arises. A bound loose by a factor of 1.7 is a working number: allocate it, run, and use 60% of it.

The difference is not that one theorem is better than the other. It is that the growth bound is over all matrices and the structural bound is over all permutations of one matrix, and the second set is enormously smaller. A bound is as tight as the set it quantifies over is small, which is obvious in retrospect and is a useful thing to ask about any bound before deciding what to do with it.

It also suggests the reading that makes both of them make sense at once. The growth bound says elimination is safe; the structural bound says the factor will fit. Neither says what will happen. For that there is only the measurement, which is why every figure in this field prints one.