Numerical Linear Algebra
where the algebra stops being the arithmetic
Theme

The thread: Where the algebra stops being the arithmetic — page 8

Essays 169 to 192 of the 271 on this theme, in the same order.
5860626466687010⁻¹⁷10⁻¹⁶terms in the dot productmean relative errorthe kernel changes herea constant in a libraryone accumulator3.2·10⁻¹⁷four accumulators2.1·10⁻¹⁷step at the cutoff1.6cutoff64the problem did not changethe loop did The answer that depends on the machine

The length that changes the kernel

A dot product's accuracy steps by a factor of 1.57 between 63 and 64 terms, on vectors drawn identically at both lengths. Nothing about the problem changes there. A library switches from one accumulator to four, at a constant in somebody else's source file.

 05101520253035059118177236295vertices eliminatededges of fill so farminDegree: 71natural: 125reverse: 125random: 160maxDegree: 293fill, by orderingminDegree71natural125reverse125random160maxDegree293edges to start60eliminating a vertex makes a cliqueand the order decides how big The matrix that is a graph

Eliminating a vertex is a graph operation

Gaussian elimination on a Laplacian deletes a vertex and joins its neighbours into a clique with conductances wᵢwⱼ over Σw. The matrix that remains is still a graph — symmetric, zero row sums, nonpositive off the diagonal — and the ordering decides whether the fill is thirty-one edges or four hundred and sixty-five.

 12345678910111210⁻¹⁵10⁻¹³10⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹degree of the approximantworst relative error on the target setchosen by the residualdegree reached12error there10⁻¹²even support, same degree1.1·10⁻⁷advantage702nearest support point0.096clustering ratio1.5filled: points chosen by the erroropen: points spread evenly The eigenvalue problem that is not linear

The points the algorithm chose

A rational approximant whose support points are picked by its own residual clusters geometrically at a branch point nobody named — recovering by measurement the rule a hand-built approximant is given. At degree ten it is seven hundred and fifty times more accurate than the same form with its points spread evenly.

 57911131517192110⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²degreerelative error of the worst rootmeasuredκ × ua polynomial is its coefficientsworst root, measured0.0076predicted, κ × u0.0023root condition number10¹³largest coefficient1.4·10¹⁹the roots are integersand the coefficients are not the roots Two errors, and whose fault they are

The roots are not the coefficients

A polynomial whose roots are the integers one to twenty, expanded exactly, handed to the routine every library uses. The computed roots are wrong in the third digit, the computation is backward stable for the matrix it factorised, and above degree eighteen the coefficients are not double-precision numbers at all.

 12345678110¹partial orders kept at each pairingexcess over the exhaustive orderthe exhaustive orderdashes above: each rule's worst casepair by cheapest productpair by smallest resultone knob, two rulesby cost, width 11.5by cost, width 21.1by size, width 11both, width 81by size, worse wider24by cost, worse wider3the width repairs the bad ruleand does not improve the good one Where the flop count stopped predicting the time

The search that got worse as it widened

Keeping two candidate orders instead of one removes a third of the cheapest-product rule's excess, and by eight it has removed all of it — the two greedy rules this field separated become the same rule. On twenty-four of sixty networks a wider search returns a dearer order than a narrower one, and it is the better rule that it more often makes worse.

 10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²1λ, the ridge on every subproblemerror, and the largest rank-one termthe terms a swamp was growingthe error it costs a fit that has an answerthe error it costs one that has nota bound, at its own priceterms, no ridge6.7terms, heaviest1.2sweeps, no ridge3000sweeps, heaviest69benign error at λ0.089terms still the right terms1the plateau is boundedand the answer is biased by exactly λ When the index is a tuple

The repair that costs exactly itself

A ridge on every subproblem is the standard cure for a swamp and it works — the terms stop at 1.61 instead of climbing past 6.7, and a run that never finished finishes in 798 sweeps. On a tensor that does have an answer the error it costs is the ridge itself, to within a factor of two, at every setting from 10⁻¹⁰ to 10⁻¹. And the terms it returns are still the right terms.

 what each reachesbest pair0.0034‖x‖ alone0.0048‖L₁x‖ alone0.0035λ₁, on ‖x‖λ₂, on ‖L₁x‖ — log₁₀0-6-5-4-3-2-100-6-5-4-3-2-10outlined row and column: one penalty switched offred cell: the best pair, 0.00343darker is worse, on a logarithmic scalethe optimum sits on or beside an edge Regularisation, and the answer that is chosen

A second penalty is not a second parameter

Penalise ‖x‖ and ‖L₁x‖ at once and there are two λ to choose. Over ten draws on five signals the best pair beats the better single penalty by between 0.00% and 3.1%, and one of its two parameters is exactly zero on 30 to 70 per cent of draws. Choosing the wrong one of the two costs up to 54%. The surface is a choice between two curves with a knob nobody needs.

 an ordering that is theregreedy at n = 764best ordering at n = 72its largest multiplier1456710¹10²matrix size ngrowth factorgreedy, ties forwardsgreedy, ties backwards, ε = 10⁻¹²greedy, ties backwards, ε = 0the best orderingthe best order is 2 3 4 5 6 1 at n = 6one cyclic shift of the rows Elimination, and the swap

The order the greedy rule cannot choose

Wilkinson's matrix is the standard demonstration that partial pivoting's growth bound of 2^(n−1) is attained. It is attained by the row order the greedy rule picks, and not by the matrix: a single cyclic shift of the rows gives growth 2 at every size, with no multiplier above one. At n = 7 that is 64 against 2. And perturbing one entry by 10⁻¹² leaves the good order exactly where it was while putting every tie-break of the greedy rule back on 64.

 what a block storesnumbers in T, block of 11numbers in T, block of 16136‖QᵀQ − I‖ there3.9·10⁻¹⁵10⁻¹⁶10⁻¹⁴10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶block size‖QᵀQ − I‖124816T perturbed by 10⁻⁶T as computed1, 3, 10, 36, 136 computed numbers in the triangleand a factor of five in what they produce Orthogonality, measured

A triangle where the scalar was

Every level-3 QR assembles a block of reflectors into Q = I − Y T Yᵀ, and T is computed by a recurrence whose inputs are its own previous columns. A block of sixteen carries 136 computed numbers where sixteen separate reflections carry sixteen. The orthogonality it produces is 3.9·10⁻¹⁵ against the single reflector's 7.8·10⁻¹⁶ — a factor of five for a hundred and thirty-six times as many things that have to be right.

 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 Sparsity, and what elimination costs

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.

 six tolerancesblocks, every tolerance250slope per unit of rank5.1worst residual3.6·10⁻¹⁵234567010203040mean rank of the low-rank blocksnumbers per doubling2 digits4 digits6 digits8 digits10 digits12 digitsthe partition is the same at every toleranceand the rank rises one per two decades Neither sparse nor dense

The partition that does not move

The storage curve's growth rate more than triples between two digits and twelve, and the earlier measurement said so without saying what moves it. There are two candidates and the measurement separates them decisively: the partition is identical at every tolerance — 250 blocks, 156 of them low-rank, 94 dense — and the mean rank rises by exactly one per two decades of accuracy. The slope is 5.11 numbers per unknown per doubling per unit of rank, with a worst residual of 0.03 over six tolerances.

 the problem you posedA = H10b = A·(1, 2, …, 10)the problem it answered exactlyA + δA, b + δb‖δ‖ / ‖A‖ = 2.3·10⁻¹⁷the answer you wantedx = (1, 2, …, 10), exactlythe answer you gotx̂, wrong by 2.7·10⁻⁴ relativebackward error 2.3·10⁻¹⁷forward error 2.7·10⁻⁴κ = 1.6·10¹³κ · η = 3.6·10⁻⁴, and the measured forward error is 2.7·10⁻⁴.The algorithm is not at fault. The problem is.H10, LU with partial pivotingresidual and error differ Two errors, and whose fault they are

Three errors and one number

This site's identity has two factors and a division of blame between them. Two fields have now added a third party and a fourth, and only one of the four is a property of anything — the others are decisions, made before the arithmetic, reported by nothing.

 10⁻²10⁻¹110¹10⁻¹⁷10⁻¹⁵10⁻¹³10⁻¹¹|λ|backward error of the computed pairone residual, four denominatorspairs measured16worst, all three1.4·10⁻¹⁵worst, K alone1.8·10⁻¹⁴worst, M alone9·10⁻¹²K-only factor, low1.1high46solid: every coefficient may movedashed: only one may The eigenvalue problem that is not linear

A perturbation that moves every coefficient

The backward error of a polynomial eigenpair is measured against perturbations of all three coefficients at once. Restrict it to the one coefficient anybody is willing to move and the same computed answers are stable at one eigenvalue and unstable at another, by a factor that runs from 1.06 to 6,370 across a single spectrum.

 10⁻¹110¹10²10³1234eigenvalue of the preconditioned Laplacian, and the stretchbreadth-firstκ = 34.62randomκ = 55.17random, secondκ = 37.11the bar is the spectrumthe tick is a combinatorial bound The matrix that is a graph

A preconditioner that is a tree

Every eigenvalue of a tree-preconditioned Laplacian is at least one and at most the total stretch — a combinatorial integer with no arithmetic in it. Measured, the bound is two to four times loose, and on a grid the preconditioner makes the conditioning worse by a factor of 1.85 at every size.

 81930415263110²10⁴10⁶10⁸10¹⁰10¹²size of the matrix|rₙₙ| ÷ σₘᵢₙthe two agreec = 0.2c = 0.35c = 0.5no ceilingratio at n = 1021ratio at n = 646.8·10¹⁰interchanges, anywhere0the greedy rule never had a choiceand the gap grows with every row Eigenvalues, singular values, rank

The cheap rank and what it cannot see

Almost nobody computes singular values to decide a rank. The standard substitute is QR with column pivoting, read off the diagonal of R — and there is a triangular matrix on which the greedy rule makes no interchange at all, has no better column available at any step, and reports a matrix eight orders of magnitude further from singular than it is.

 01210¹10²10³the three kinds of runsweep at which the test firesno answer to reachbadly conditionedordinaryhollow: never fired, drawn where the run endedwhat it fires on, and whenthreshold0.25boundary, fired6boundary, median sweep22collinear, fired6collinear, median sweep759ordinary, fired0the value does not separate themand the sweep does When the index is a tuple

The test that is a deadline

The quantity that separates a boundary from slow convergence fires on every run that has nothing to reach, at sweep 21 of four thousand, and on no ordinary fit. It also fires on every badly conditioned fit that does converge — at sweep 758. No threshold between 0.10 and 0.40 separates those two, and the sweep it fires at separates them by a factor of thirty-six.

 n = 512, ε = 10⁻⁸least storage, at leaf16numbers per unknown there133at the largest leaf211050100150200leaf sizestored per unknown48163264stored per unknownshare in the dense blockslarge dot: the least storage availableand the dashed curve is why there is one Neither sparse nor dense

Two knobs on one number

The tolerance raises every block's rank and leaves the partition alone. The leaf size does the opposite — it changes how many blocks there are and does not move a single rank. Both act on the same storage, and only one of them has a best setting: at eight digits the least storage is at a leaf of 16 and the largest leaf tried costs 59 per cent more. The best leaf rises with the accuracy, from 4 at two digits to 16 at twelve, and the influence runs only that way.

 10⁻⁹10⁻⁸10⁻⁷10⁻⁶10⁻⁵10⁻⁴10⁻³10⁻²110¹10²10³10⁴standard deviation of the noise σmedian growth factorthe marginevery entrynonzeros onlyrelative, per entryone error in Estep 0.3, margin 0.00564, thirty drawsno noise1.1·10⁴every entry, σ = 10⁻⁵536nonzeros only, σ = 10⁻³1.1·10⁴one error in E, σ = 10⁻², draws kept21every curve is the same matrixonly what the noise touches changes Elimination, and the swap

Noise the growth amplifies

The shooting matrix's growth of 1.1·10⁴ fell under noise as its pivot margin divided by the noise, and the explanation offered was noise reversing partial pivoting's choices. But noise of 10⁻⁵ is five hundred times smaller than that margin. Add the same noise only to the entries that are not zero and every draw keeps the whole growth up to 10⁻³. The dense noise was not reversing the comparisons by itself: it sat in the zeros, the elimination multiplied it by the growth already made, and a comparison flips when that product reaches about four tenths of the margin — at every noise level from 10⁻⁶ to 10⁻³.

 10⁻¹110¹10²s, the smallest interpolation pointrightmost pole of the reduced model−10−10+1+10unstable above zeronumerical range edge, -0.986system's abscissa, -3.49two-sided, order 3one-sided, order 3one-sided, order 6Pe 10, nineteen placementssystem's rightmost eigenvalue-3.5edge of the numerical range-0.99two-sided, order 3: unstable placements4one-sided, order 3: unstable placements0one-sided, order 6: unstable placements0signed logarithmic axisshaded: poles in the right half plane Reduction, and what a model is for

Half the conditions and a certificate

A two-sided interpolatory reduction of a stable convection–diffusion system came back unstable at four placements of nineteen. The one-sided reduction built from the same kind of solves came back stable at all 152 placements measured across eight Péclet numbers, and not by luck: every pole of a Galerkin model lies inside the numerical range of the operator, whose edge here is exactly the diffusion term's largest eigenvalue, −0.986 at Péclet 10. The price is the slope conditions, and spent as six one-sided shifts instead of three two-sided ones it beats the stable two-sided model at seven placements of ten.

 1220283644526068768492100101214161820grid points nrelative error against the continuous signal, per centsampled kernel, hatsintegrated, hatsintegrated, splinespline, no λ1% noise per sample, median of sixteen drawssampled kernel, hats: 16 points0.18integrated, hats: 16 points0.16integrated, spline: 16 points0.16sampled kernel, hats: 96 points0.14integrated, hats: 96 points0.14integrated, spline: 96 points0.14spline, no λ, its best grid0.15dotted: the 96-point answerevery point is a grid with its own λ Regularisation, and the answer that is chosen

The grid on which the discretisation stops mattering

Without regularisation, integrating the blur's kernel against cubic splines beat sampling it at 0.1% noise and lost to it at 1%. Give each discretisation its own best λ on every grid and the difference shrinks to nothing where grids are fine — 0.17, 0.28 and 0.20 per cent apart on 96 points at the three noise levels, with every discretisation choosing the same λ — and stays at 17 to 18 per cent on 16 points. The choice between them is a choice of how coarse a grid can be: at 0.1% noise the integrated discretisations reach the fine-grid answer on 26 points and the sampled one needs 40.

 00.250.50.75100.250.50.751xuexactcentral differencesupwindthe oscillation is exact‖Ax − b‖/‖b‖ for the central answer4.6·10⁻¹⁸values outside [0, 1]16worst excursion0.52the dashed lines are 0 and 1, which the equation guaranteesno solver was involved Iterating, instead of factorising

The stencil that is not symmetric

Past a cell Péclet number of exactly one — measured by bisection at 1.0000000000000002 — the central-difference solution of a convection–diffusion problem oscillates from point to point and leaves the interval the equation guarantees, at 16 of 31 grid points. It is the exact solution of its own linear system, to 4.6·10⁻¹⁸. No solver was involved.

 110¹-101pieces the inner product was summed inqᵢᵀqⱼ ÷ 2.78·10⁻¹⁷the true valuethe norm and its partsdistinct values5of runs10pairs with no fixed sign125of pairs1128κ of this inner product1.1·10¹⁷‖QᵀQ − I‖ moves by1the aggregate is stableand no entry of it is The answer that depends on the machine

An inner product with no fixed sign

‖QᵀQ − I‖ is how this site turns "orthogonal" into a number, and across ten partitionings it moves by 2.4%. The entries it is built from are not so lucky: 125 of the 1,128 off-diagonal pairs take both signs, and one of them takes five different values including zero.

 57911131517192110⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹vertices of the complete graphrelative sizethe algorithm is not the problemsizes counted18worst relative error9.4·10⁻¹⁵first wrong at n =16digits in the answer17spacing of doubles there8largest integer a double holds9·10¹⁵the elimination stays at the rounding levelthe format runs out The matrix that is a graph

A count that comes out of a determinant

The number of spanning trees of a graph is the determinant of its grounded Laplacian, so it is a whole number known in advance. The elimination that computes it is backward stable at every size — and from sixteen vertices the answer is wrong, because the count has seventeen digits and a binary64 has sixteen.

 22.539343.078693.618034.157384.6967210⁻¹110¹10²stiffness damping βoverdamping marginβ* = 3.23607two routes to a boundaryclosed form β*3.2by certificate3.2difference4.7·10⁻¹³bisection steps44a factorisation that completesand a sine, agreeing to twelve digits The eigenvalue problem that is not linear

A class a longer chain takes away

Symmetry survives a bigger problem. Hyperbolicity does not. The damping that certifies a chain of seven masses is refused by a chain of eight, the damping the class demands grows like the length without bound, and the certificate has to be earned again at every size — which costs one Cholesky, and the alternative is a proof quietly inherited from a smaller problem.

All themes