Numerical Linear Algebra
where the algebra stops being the arithmetic
Concept

Sparsity — where it appears

The property of having mostly zero entries, which decides what a factorisation costs and disappears from every factor unless the ordering protects it. It decides what a factorisation costs and it does not survive the factorisation, which is why the ordering that protects it is a decision nothing numerical makes.

Named by 20 essays across 7 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 051015202530354010⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹singular value, in orderσ ⁄ σ₁eight digitsan independent draw per entry1 ⁄ ra cliff, and a control1/r rank at 10⁻⁸5log r rank at 10⁻⁸5noise rank at 10⁻⁸96σ₂ ⁄ σ₁0.024σ₆ ⁄ σ₁2.8·10⁻⁹the block has full rankand five useful columns

A block nobody can call sparse

A 96 × 96 block of a kernel matrix has ninety-six nonzero singular values and five that matter. It has no zero entries, it is not described by fewer numbers than it contains, and neither of the two ways this collection already knows to make a large matrix affordable applies to it.

hierarchy · Off-diagonal rank 0408012016020024010⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹iteration‖e‖ ⁄ ‖e₀‖ in the A-normmeasuredκ bound119 steps40×40, spectrum spread evenly in logbound permits 1417

The rate the condition number predicts

Conjugate gradients converge at a rate governed by the square root of the condition number. That is a bound rather than an estimate, it is provable, and it is loose enough that provisioning iterations from it wastes nine out of ten.

iterative · Krylov natural1739reverse Cuthill–McKee1354minimum degree1026nested dissection1413matrix: 408 entries · dense factor: 10440bandwidth 12 · 4.26× the matrixbandwidth 12 · 3.32× the matrixbandwidth 123 · 2.51× the matrixbandwidth 108 · 3.46× the matrixn = 144, five-point stencilevery ordering fills in; none avoids it

The order decides the memory

Four elimination orderings on one matrix give factors of 1,739, 1,354, 1,413 and 1,026 entries. All four factorisations are exact, all four return the same answer, and the one with the better asymptotics is not the one that wins.

sparsity · Ordering the matrix43 entriestip eliminated first253 entriestip eliminated last43 entries‖A − LLᵀ‖/‖A‖, tip first1.4·10⁻¹⁶‖A − LLᵀ‖/‖A‖, tip last0dense factor is n(n+1)/2 = 253 · sparse factor is 2n − 1 = 43one row swapped to the endnothing numerical chose between them

Two ends of the same arrow

One matrix, one row moved from the front of the elimination order to the back, and the factor goes from completely dense to no fill at all. Both factorisations are exact to rounding, and nothing numerical chose between them.

sparsity · Fill 051015202530354010⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹iteration‖r‖ / ‖b‖plain CGIC(0) CGwhat the preconditioner didκ(A)48κ(L⁻¹AL⁻ᵀ)5.1‖A − LLᵀ‖/‖A‖0.0832D Laplacian, n = 100√κ ratio predicts 3.07×

Changing the condition number on purpose

Preconditioning is usually introduced as a trick that makes an iteration converge faster. It is not a trick. It is solving a different system with the same solution and a condition number chosen rather than inherited, and the new condition number is computable.

iterative · Preconditioning 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 10⁻³10⁻²10⁻¹1110¹10²10³10⁴pivot threshold τgrowth factor · entries in L+U, ÷ entries in Agrowthfillthe library default‖PA − LU‖/‖A‖ at τ = 0.12.9·10⁻¹⁶growth at τ = 0.138entries in L+U at τ = 0.1372one knob, two measurements, opposite directionsand the default is most of both

A threshold between fill and growth

One number decides how small a pivot an elimination will accept. At 0.001 the factor holds 172 entries and the matrix grows by 1,330; at 1 it holds 260 and grows by 1.2. The libraries ship 0.1, and the measurement says why.

sparsity · Sparse pivoting 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 κ(A) = 10⁶ throughout · κ(H) = 100 · the answer is the same answer for every basisorthonormal — κ(Z)1κ(ZᵀHZ)25.6relative error1.07·10⁻¹⁵first m basic — κ(Z)1.99·10⁸κ(ZᵀHZ)3.8·10¹⁶relative error0.0518pivoted basic — κ(Z)2.06κ(ZᵀHZ)31.9relative error6.71·10⁻¹⁶what the choice costsdensity, orthonormal1density, fundamental0.5κ(ZᵀHZ) ÷ κ(Z)², naive0.96error, pivoted choice6.7·10⁻¹⁶every one of them is a basisand one of them loses fourteen digits

The basis nobody chose on purpose

A method that eliminates a constraint has to pick a basis for its null space, and every basis is correct. Their condition numbers are eight orders apart, the reduced problem inherits the square, and the choice is usually made by a one-line rule nobody thought of as a numerical decision.

orthogonality · Null-space basis 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 A6×6B6×6I ⊗ A + Bᵀ ⊗ I36×36one equation, two objectsentries in A and B72entries in the coefficient matrix1296two routes, relative gap1.7·10⁻¹⁶‖AX + XB − C‖/‖C‖1.4·10⁻¹⁶the small squares are the problemand the large one is the notation

The elimination the matrix does not need

The Kronecker form of AX + XB = C is dismissed with a hundred million entries and (2/3)n⁶ operations. Both price an elimination, and after the reduction both routes take, the matrix has exactly n³ nonzeros, none of them above the block diagonal, and nothing left to eliminate.

structure · Matrix equation 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 1357911131510⁻¹⁶10⁻¹⁴10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²1k, where h = 1 − 10⁻ᵏrelative error of 1 − hsubtractwith a correction stepcorrection kept as a pairfrom the reflectors40 × 6k = 12, subtracting1.1·10⁻⁴k = 12, with a correction1.1·10⁻⁴k = 12, from the reflectors3·10⁻¹⁵reflector operations900dashed: a unit of roundoff over the divisorthree routes sit on it and one does not

The factor a sparse code keeps anyway

Every deletion diagnostic divides by one minus a leverage, and computing it as a subtraction loses a digit for every decade the leverage is from one. The route that does not subtract needs the orthogonal factor, which a sparse factorisation is supposed not to have. Three repairs that avoid it all fail at exactly a unit of roundoff over the divisor — and the fourth, which reaches the orthogonal factor through the Householder vectors a sparse code keeps in order to solve anything at all, returns the same bits as a stored factor in 900 operations.

leastsquares · Leverage a 28 × 28 saddle-point matrixnatural order, entries113sparsest order, entries72two-by-two pivots, loosest test0and at the strictest3023466992115the constant in the pivot testentries in the factor0.10.30.640.80.95133natural ordersparsest ordergrowth, naturalgrowth, sparsestnumbers above the points: two-by-two pivots takenboth curves favour the sparser order

The freedom a symmetric factorisation does not have

Permuting rows and columns together leaves no column to choose, so the conflict between the sparsest pivot and the sound one should be worse rather than better. On a saddle-point matrix whose constraint rows have no diagonal entry at all, it is not there: taking the sparsest available pivot holds 70 entries against the natural order's 113 and a growth of 1.28 against 1.83 — better on both currencies at once, at every setting of the pivot test. The two-by-two blocks that make it legal cost 1.33 entries apiece.

sparsity · Sparse pivoting 0481216202410⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹matrix–vector products, mrelative error in eᴬbforming e^A: 3·10⁻¹⁶crosses at m = 18the vector, not the matrixsteps to the dense answer18dimension100Krylov megaflops0.36dense megaflops2the exponential that is computedis 18×18

The vector was what was wanted

Nobody who computes a matrix exponential wants the matrix. They want eᴬᵗb — one vector, the state of a system at a later time. Twenty matrix–vector products get it to sixteen digits on a hundred-by-hundred problem, without ever forming a hundred-by-hundred exponential, and the exponential that does get computed is twenty by twenty.

spectra · Matrix function 10⁻¹⁷10⁻¹⁵10⁻¹³10⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹10⁻¹⁴10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²1εrelative error in J(x)vforwardcentralcancellationtruncationagainst a derivative that is exactforward floor1.3·10⁻¹⁰central floor1.1·10⁻¹²truncation slope, forward1truncation slope, central2no ε reaches the roundoffand the analytic derivative is free of the choice

An operator with no entries

At the sizes where linear algebra is expensive the matrix does not exist. What exists is a subroutine that returns Av. Every Krylov method survives that unchanged; every algorithm that reads an entry disappears. And the derivative such a code computes is accurate to ten digits instead of sixteen, which turns out to cost nothing at all.

iterative · Matrix-free

Named alongside it

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

Fill-inFill-reducing orderingSymbolic factorisationGaussian eliminationPermutationCondition numberFlop countMinimum degreeResidualElimination graphGrowth factorNested dissection

All concepts