Numerical Linear Algebra
where the algebra stops being the arithmetic
Concept

Symbolic factorisation — where it appears

The phase that computes a sparse factor's nonzero pattern from the matrix's graph before any arithmetic. What it produces is a bound rather than a prediction, since an entry it allocates may turn out to be numerically zero.

Named by 13 essays across 4 fields — each of them below, with the objects they name alongside it.

the matrix, lower triangle408 entriesits Cholesky factor1739 entries · 1331 created→‖A − LLᵀ‖/‖A‖1.4·10⁻¹⁶fill, symbolic1331fill, numeric1331n = 144 · density 3.2% · bandwidth 12same matrix, renumberedthe answer is identical to rounding

The factor is not sparse

A sparse matrix has a factor that is not sparse, and the gap between them is the entire reason iterative methods exist. The entries elimination creates can be counted before any arithmetic runs, from the graph alone.

sparsity · Fill the matrix105 entriescorner first — sparsest227 entries, growth 1.9·10¹¹largest first — safe242 entries, growth 1.19the middle factor is the smaller one, and its answer has no correct digitsboth factorisations reproduce the matrix‖PA − LU‖/‖A‖, sparsest3.8·10⁻¹⁷‖PA − LU‖/‖A‖, pivoted5.4·10⁻¹⁷forward error, sparsest3·10⁻⁵forward error, pivoted4.8·10⁻¹⁶red marks are entries elimination createdthe fill argument and the stability argument disagree

Structure and stability stop being separable

The sparsest variable to eliminate on this matrix has a diagonal entry of 10⁻¹². Eliminating it produces the smaller factor, reproduces the matrix to 3.8·10⁻¹⁷ — better than pivoting does — and returns an answer wrong in the fifth digit.

sparsity · Sparse pivoting 545545545455554545545545545454555555454545545545545455554545545545rows and columns, in the order the points arrivedark: kept dense · light: two thin factors, rank printedthe strong partitionblocks112kept dense46largest rank5numbers stored2.7·10⁴relative compression error3.4·10⁻¹⁰the picture is decidedbefore a number is read

Which pairs are allowed to be small

A hierarchical representation is a partition of the matrix into blocks, and the rule that produces it reads four numbers per pair of index clusters and not one entry of the matrix. On a 256-square it yields 112 blocks, 66 of them stored as two thin factors, none of rank above five.

hierarchy · Admissibility the bound211no pivoting127τ = 0.00198τ = 0.003100τ = 0.01106τ = 0.03111τ = 0.1116τ = 0.3117τ = 1138entries in Ua bound, and its slackthe bound, from the graph alone211the worst that occurs138loose by1.5no arithmetic was done to compute the bound — only the pattern of AᵀA and its elimination graphGeorge and Ng: U fits inside chol(AᵀA), whatever the swapsprovable, cheap, and loose

What the symbolic phase can only bound

Without pivoting, the fill can be computed from the graph and the count is exact — 233 predicted, 233 measured. With pivoting it is 233 predicted and 242 measured, and what survives is a bound that is right at every threshold and loose by 1.7 times at the largest grid drawn.

sparsity · Sparse pivoting -18-15-12-9-6-300log₁₀ ‖PKPᵀ − LDLᵀ‖ / ‖K‖share of orderingsbestworstexistence and stabilityfactorise1unregularised0.69worst growth6.4·10⁵growth × γ0.64the ordering is free to chooseand not free of consequence

The regularisation that legalises every order

Perturb a saddle-point matrix's two blocks in opposite directions and it acquires a factorisation with a diagonal D under every symmetric permutation — not under a good one, under all of them. Five hundred random orderings, five hundred successes, and a growth factor that spans six orders across them.

constraint · Quasi-definite 02468101210⁻¹⁸10⁻¹⁶10⁻¹⁴10⁻¹²10⁻¹⁰10⁻⁸member of the sequencebackward error of the solvefresh order · equilibrated kept orderringed: another pivot replaced hereone order, two sets of unitsfresh order, last member9.5·10⁻¹⁷kept order, last member4.8·10⁻⁹equilibrated, last member6.5·10⁻¹⁷pivots replaced, kept2pivots replaced, equilibrated0the perturbation floor √u1.1·10⁻⁸a reused pivot order is safe or notdepending on what the rows are measured in

The order that was right last time

A pivot order computed once and reused across a sequence saves the symbolic phase, and the price is that a pivot which was large may now be small. Replacing it with √u·‖A‖ costs eight orders of backward error and iterative refinement recovers a factor of 8.8 of them. Divide each row by its largest entry first and the same reuse costs nothing at all.

sparsity · Sparse pivoting -1-0.75-0.5-0.25010⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹10⁴log₁₀ σ — the barrier's reduction factorresidual after one reused stepconvergedthe pattern free, the factors notentries moved6off-diagonal0survived at σ = 0.996survived at σ = 0.106 steps0 stepsthe few entries that movedare the ones that dominate

What survives one step of the barrier

An interior-point method solves the same system dozens of times with the same pattern and different numbers, and exactly p entries change between one step and the next. The pattern is reusable for ever. The factorisation is reusable for none of them, and the threshold that says so is a reduction factor of about a per cent against schedules that use ten.

sequence · Reuse natural113 predicted · 113 countedminimum-degree63 predicted · 63 countedreverse Cuthill–McKee63 predicted · 63 countedthe shaded entries are fill: zeros of K that the factorisation makes nonzeroallocated before the numbersnatural113minimum-degree63reverse-cuthill-mckee63predicted minus counted0the symbolic phase decides the memoryand nothing later is allowed to argue

An ordering that does not wait for the numbers

A sparse factorisation's memory is decided by an ordering computed from the graph, and its stability by pivots computed from the values, and the two decisions fight. On one family of matrices they do not — the ordering can be chosen for fill alone, and the fill the symbolic phase predicts is the fill the factorisation produces — exactly, not as a bound.

sparsity · Sparse pivoting 6101418222610³10⁴10⁵grid side karithmeticND, totalMD, totalMD, critical pathND, critical pathdashed: the time on unbounded processorsthe same time, bought with more work

An ordering that buys processors, not time

Nested dissection loses to minimum degree on fill and on total work at every grid either measurement could draw. Read along the elimination tree a parallel factorisation works on, it does not win back the time either: its critical path is within 28 per cent of minimum degree's at every size from 8 to 24 points a side, in both directions, and the tree heights and widest columns are nearly the same. What it wins is the ratio. Its total work over its critical path — the most a factorisation on unbounded processors can speed up by — grows from 2.7 to 4.5 while minimum degree's stays between 2.1 and 2.6.

sparsity · Ordering least possible76minimum degree76worst of 300 tie-breaks81reverse Cuthill–McKee85natural order9920 unknownstie-breaks at the least0.44minimum degree ÷ least1every order searched, by the set already eliminatedthe greedy rule found the minimum

The least fill there is

Finding the elimination order with the least fill is NP-hard, and that is a statement about the hardest graph and the largest size. On a graph of twenty vertices every one of the 20! orders can be searched at once, through the million sets of vertices already eliminated, and the least fill is a number. On the 4×4, 4×5 and 3×7 grids minimum degree finds it exactly. On eighty random sparse graphs of eighteen vertices it finds it on 53 and misses by at most 7.6 per cent, and on every one of the eighty some breaking of its ties finds it.

sparsity · Ordering all 2^16 subsets searchedboth minima, one order0.97minimum degree optimal0.82worst fill overshoot1.1worst work overshoot1.111.021.041.0611.041.081.121.16fill ÷ the least filloperations ÷ the leastovershoot squaredon the diagonal: the two overshoots are equalon the upper curve: the arithmetic overshoot is the square

Two minima that are one minimum

The order that decides the memory found the operation count behaving like the square of the fill, which leaves room for an order with slightly more fill but a shorter heaviest column to do less arithmetic. Searched exactly over every elimination order on forty graphs, that order does not exist: one order attains both minima on thirty-nine of forty, and on the fortieth the least-fill order's arithmetic is 1.0099 times the least. Minimum degree attains both on the same thirty-three graphs and neither on the same seven.

sparsity · Ordering total operationswork, no dissection10⁵work, full dissection1.5·10⁵work, one level1.3·10⁵on unbounded processorscritical path, no dissection4.1·10⁴least, at depth 43.4·10⁴critical path, one level5.4·10⁴10⁵switching depthoperations01234568total operationscritical pathdepth 0 is minimum degree, depth 8 is nested dissectiondepth 1 is worse than both

The depth that is worse than both ends

Nested dissection to a chosen depth and minimum degree below it is the ordering codes ship, and sweeping the depth was supposed to find a setting that keeps most of dissection's parallelism for most of minimum degree's work. It does not exist: total work rises with the depth at every grid size, and one depth — the first — is worse than both extremes on work and on the critical path at all four sizes measured. One bisection buys nothing because there is no recursion under it to amortise the separator.

sparsity · Ordering 80 matrices, τ = 0.1, ties to the largestno column search, median fill244one column136every column109share of the gain at one column0.804896144192240columns the search may look atmedian fill123468121624allno column search at allmedian fillmedian growthworst growth — see the captionthe fill is the quantity the width buysthe worst case is the draw

How few columns the search needs

A full row-and-column pivot search is quadratic in the active submatrix at every step, and no library performs one. Looking at a single sparsest column takes the median fill from 244 to 136 where the full search reaches 109 — four fifths of the benefit for a linear scan — and that share is 79, 83, 80, 89 and 92 per cent across five thresholds. The worst growth appears to favour the narrow search by a factor of six, and on the next draw it favours the wide one by two.

sparsity · Sparse pivoting

Named alongside it

The objects these essays reach for when they reach for this one.

Fill-inSparsityFill-reducing orderingGrowth factorMinimum degreeElimination graphElimination treePermutationFlop countGaussian eliminationIterative refinementNested dissection

All concepts