The collection

Every essay — page 14

Essays 313 to 331 of 331, in the same order.

The answer that depends on the machine

Every other field here asks how wrong an answer is. This one asks how many answers there are. A parallel reduction adds a vector up in however many pieces there are workers, in whatever order they finish — so the same program, on the same data, at the same precision, returns a different number on a different machine, and every one of those numbers satisfies the published bound. The disagreement is a quarter of κu and the bound is ten thousand times larger, which is why nothing reports it. It matters where a number is compared to something: a stopping test, a rank test and a definiteness test each turn a real number into a verdict, and a verdict has no last digits for a disagreement to hide in. One matrix here has three different numerical ranks and one solve has thirteen different bills. And the smallest instance needs no parallelism at all — a multiply the compiler was allowed to fuse, which is one rounding, and which decides the sign of a determinant whose value is one.

-12-10-8-610⁻¹³10⁻¹¹10⁻⁹10⁻⁷log₁₀ of the residual tolerance asked forforward error of the answer1.34×1.48×1.71×1.17×the band does not closetolerances swept4runs at each7accuracy gained1.5·10⁶ratio at 10⁻⁶1.3ratio at 10⁻¹²1.2the bars falland they keep their height

The tolerance that buys no agreement

Ask for four more orders of accuracy and you get them — the answers improve by a factor of 1.5 million. The ratio between the best and the worst run is 1.34, 1.48, 1.71 and 1.17 across the same sweep. The band falls and it does not close.

4 figures · Stopping test rung 2
13579111310⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹index ksingular valueσ₁ · 10⁻¹⁴, the thresholdone matrix, three rankspartitionings7lowest rank10highest rank12threshold7.3·10⁻¹⁴σ₁7.3the curves separate in the noiseand the threshold is drawn through it

A rank that depends on the thread count

One 60 × 14 matrix, one threshold, seven partitionings of the inner products that build its Gram matrix — and numerical ranks of 12, 12, 12, 10, 10, 11 and 11. Not a digit of an answer: the number of columns a model built from this matrix would have.

5 figures · Rank rung 4
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 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.

5 figures · Algorithm selection rung 1
-15-13-11-910⁻¹³10⁻¹¹10⁻⁹10⁻⁷log₁₀ of the defect, relative to the entry it sits inhow far the answer movedthe machine, on its owna tolerance with two sidesmachine band3.2·10⁻¹²smallest caught10⁻¹²its nearest run4.6·10⁻¹²defects hidden3window, factor1.4below the band nothing is visibleand above it everything is

What a regression test can ask for

The machine's own variation on one solve is 3.2·10⁻¹², and the smallest defect whose answers clear it is one part in 10¹². The tolerance exists, it is bracketed on both sides by a factor of 1.42, and it is neither zero nor the 10⁻⁸ that usually gets typed.

4 figures · Regression tolerance rung 1
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

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.

5 figures · Orthogonality rung 2
00.10.20.30.40.50.60.70.80.9110⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹decay rate of the spectrumrelative spread of the rank-k errorthe seedthe machinetwo kinds of variationseeds drawn8partitionings7seed spread0.31machine spread9.8·10⁻¹⁶ratio3.1·10¹⁴one of these is recordedand it is the large one

The variation that comes with a seed

A randomised low-rank approximation's error moves by 31% between draws and by 10⁻¹⁵ between partitionings of one draw. In the one field on this site whose answer already comes as a band, the machine is inside the width of the line — and it is still there.

5 figures · Run-to-run variation rung 1
κ of the suma component of b − Ax8.09·10¹⁷ad − bc, near-degenerate3.6·10¹⁶qᵢᵀqⱼ, an orthogonality check3.26·10¹⁶zᵀAz, a trace probe95.7pᵀAp, a curvature73.7rᵀr, a residual norm1measured, not assumedhighest8.1·10¹⁷lowest1above 10¹⁰3terms128sums of squares are safeand nobody decides anything from one

Two machines, one certificate

Nothing a solver returns says which of its answers you got. Four things could be reported instead — the summation condition number, the partition count, an exactly accumulated residual and a directed-rounding interval — and each costs about one pass over data the routine already has in hand.

4 figures · Regression tolerance rung 2
‖QᵀQ − I‖ of the implied Qκ = 10⁸, Cholesky0.37κ = 10⁸, sweep8.5·10⁻⁹κ = 10¹⁰, Cholesky1.3κ = 10¹⁰, sweep3·10⁻⁷κ = 3·10¹⁰, Choleskyrefusedκ = 3·10¹⁰, sweep6.7·10⁻⁶κ = 10¹², Cholesky1.5κ = 10¹², sweep1.4·10⁻⁴the safe run is the one that failsrefusals in the range1wholly non-orthogonal returns2the pivot it refused on-6.5·10⁻¹⁷one of these outcomes is safeand it is the refusal

The licence is not the boundary

Cholesky QR is licensed by κ²u ≪ 1, which reaches equality at κ = 9.5·10⁷ in double precision. At 10⁸ the factor it returns is already 0.37 away from orthogonal, and it goes on returning factors as far as 10¹³ — refusing at scattered condition numbers in between, at different ones for eight columns and for six.

6 figures · Algorithm selection rung 2
614223038465410⁻¹³10⁻¹⁰10⁻⁷10⁻⁴working significand bitsrelative sizeerrorresidualuthe gap the solver cannot seeerror ÷ residual at 24 bits1.5·10⁶error ÷ residual at 16 bits2.3·10⁷error ÷ residual at 8 bits2.4·10¹⁰every convergence test passesand the answer is wrong in proportion to u

The reading that never moves

Across thirty runs — five grids from 36 to 196 unknowns, six working precisions from 53 significand bits down to 8 — the residual conjugate gradients stops on stays between 1.10·10⁻¹³ and 9.95·10⁻¹³. Over the same thirty runs the error of the answer spans a factor of 2.39·10¹¹, and the step count more than doubles. The one number the run publishes is the only one that responds to neither axis.

6 figures · Stopping test rung 3

The matrix a constraint makes

Every difficult matrix in the other fields was difficult for a reason the arithmetic supplied: a Hilbert matrix arrives ill conditioned, a Wilkinson matrix grows under elimination, a kernel matrix is dense. Ask a problem to minimise something subject to a constraint and the matrix that results is difficult for a reason the algebra supplies. Its zero block is the second derivative of a Lagrangian with respect to its own multipliers, so no pivot order removes it and no precision changes that; its inertia is known before anything runs; and Cholesky does not fail somewhere on it, it fails at the first constraint row, on a number the problem already contained. Then the field's second surprise, which is the opposite one: an interior-point method drives the condition number of this matrix to 10¹⁵ deliberately, and the answer keeps fifteen digits — because the number that describes the error is not the one every library prints.

-1-0.582271-0.1645420.2531860.6709151.088641.506370eigenvalue4 negative10 positivecounted before it was formedpositive10negative4at zero0innermost ratio39the zero block is a theoremand so is the count either side of it

The zero that is not a missing entry

A constrained minimisation produces a matrix with a zero block, and the zero is a theorem rather than a sparsity pattern. No pivot order makes it positive definite, no precision changes that, and Cholesky does not fail somewhere on it — it fails at the first constraint row, on a number the problem already contained.

7 figures · Saddle-point systems rung 1
01234510⁻¹⁷10⁻¹³10⁻⁹10⁻⁵10⁻¹10³10⁷10¹¹log₁₀ κ(A)condition number, and relative errorκ(S)κ(ZᵀHZ)range-space errornull-space erroragainst a BigInt answerκ(S) at κ(A) = 10⁵4·10¹⁰κ(ZᵀHZ), all stops21range-space forward error5.3·10⁻⁶null-space forward error5.8·10⁻¹²both are the same algebraand only one squares

Two ways to remove a constraint

A constrained system can be reduced by eliminating the multipliers or by eliminating the constrained directions. Both give the same answer in exact arithmetic and inherit different condition numbers — one of them squares the constraint's, and the other does not contain it at all.

7 figures · Saddle-point systems rung 2
00.3670080.7340171.101031.468031.835040eigenvalue of P⁻¹Kwritten down, then computeddistinct3at 16φ computed1.6off the closed form2.9·10⁻¹⁴1 − φ1φthe preconditioner's effect is a theoremand the golden ratio is in it

Three eigenvalues, and two are the golden ratio

Precondition a saddle-point system by the block diagonal of its own two definite pieces and the preconditioned matrix has exactly three distinct eigenvalues — 1, and the two roots of λ² − λ − 1. A minimal polynomial of degree three means three steps, at every conditioning, and the preconditioner nobody can afford turns out to be the statement the affordable ones are measured against.

7 figures · Block preconditioning rung 1
012345610⁻³10⁻²10⁻¹110¹log₁₀ κ(A)eigenvalue of P⁻¹Kthe constraint is not in itnontrivial6at one8movement, six decades7.1·10⁻⁶drift at one6·10⁻⁵2m at onethe marks are a pencil that never saw Aand the lines are the preconditioned matrix

A preconditioner that need not know the constraint

Keep the constraint block exactly and replace the objective block by anything positive definite on the null space. The preconditioned matrix then has 2m eigenvalues at exactly one, and its remaining n − m are the generalised eigenvalues of a pencil in which the constraint does not appear. Sweep its condition number over six decades and they do not move in six digits.

8 figures · Block preconditioning rung 2
-14-12-10-8-6-4-2010⁻¹⁷10⁻¹³10⁻⁹10⁻⁵10⁻¹10³10⁷10¹¹10¹⁵log₁₀ μ — the barrier parametercondition number, and relative errorκ₂, condensedκ₂, augmentederror, condensederror, augmentedagainst a BigInt answerκ₂ augmented, μ = 10⁻¹⁴3·10¹⁵its relative error10⁻¹⁵κ₂ condensed2.4·10¹⁶its relative error0.31the same step, written two waysand only one of them is solvable

A condition number sent to infinity

An interior-point method manufactures an ill-conditioned matrix on every iteration, deliberately, because the separating of a diagonal is how it discovers which constraints are active. Written one way the answer keeps fifteen digits at a condition number of 3·10¹⁵. Written the other way — the way almost every code writes it — it has none left.

8 figures · Interior-point conditioning rung 1
-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.

7 figures · Quasi-definite rung 1
-5-3-11350eigenvalueHZᵀHZK4 negative — Cholesky of H stops at row 30 negative — a minimum on the constraint(10, 4, 0) = In(ZᵀHZ) + (4, 4, 0)one factorisation, no Zpositive, LDLᵀ of K10negative, LDLᵀ of K4negative in H4negative in ZᵀHZ0the count follows the reduced Hessiannot the Hessian

A minimum the Hessian cannot see

A Hessian with four negative eigenvalues can sit at a constrained minimum, and a Cholesky of it stops at the third row. One symmetric indefinite factorisation of the saddle-point matrix settles the question anyway — ten positive pivots and four negative — without a basis for the null space ever being formed. The count is exact in the algebra and blind in floating point, in a band that grows like κ(A)²; the route through the null space is blind in one that grows like κ(A).

7 figures · Saddle-point systems rung 3
-1-0.500.511.50eigenvalue of P⁻¹Ktriangular [[H, Aᵀ], [0, −Ŝ]] — 1 valuediagonal blkdiag(H, Ŝ) — 3 valuessteps to a residual of 10⁻¹⁰GMRES, triangular2MINRES, diagonal3‖P⁻¹K − I‖54computed |λ − 1| at c = 18.1·10⁻⁸one copy of each value against twoand the counts follow

One eigenvalue and two steps

Put the off-diagonal block back into a block-diagonal saddle-point preconditioner and every eigenvalue of the preconditioned matrix becomes exactly one. GMRES still needs two steps, because the matrix is the identity plus a nilpotent part of norm 54, and a computed eigenvalue at one comes back as a ring of radius 8·10⁻⁸ — the square root of the rounding, not the rounding. With an approximate Schur complement the triangular form leaves one copy of each value where the diagonal form leaves two, and the step count halves.

7 figures · Block preconditioning rung 3
0246810121410⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹interior-point iterationrelative error, and μcertified from iterate 1μthe iterate's errorthe crossover's errorone solve, checkediterations15first certified iterate1iterate error there0.22crossover error there4·10⁻¹⁴the active set arrives long before the digitsand a crossover collects them at once

The active set before the digits

An interior-point method takes fifteen iterations on a quadratic programme with forty constraints, and its iterate has eight correct digits at the eleventh. Take the constraints its diagonal calls active at the first iterate, solve the equality problem they define once, and check the answer against the conditions for optimality. It passes, to thirteen digits. The step's matrix had a condition number of 43 at that iterate, and 7·10¹⁵ at the last.

8 figures · Interior-point conditioning rung 2
-1-0.500.511.520eigenvalue of P⁻¹K1 − φ1φS never formedγ1smallest ν0.056furthest from φ, 1 − φ0.57MINRES steps11an exact Schur approximation, bought by changing Hthe golden ratio without S

Where the augmentation puts the cost

Add γAᵀA to the objective block of a saddle-point system and its Schur complement tends to I/γ, so the cheapest possible approximation becomes the right one and the golden-ratio spectrum arrives — within 7.6·10⁻⁶ at γ = 10⁶. MINRES falls from 21 steps to 6. The inner solve with the augmented block rises from 14 conjugate gradient steps to 43, their product does not fall at all, and the answer loses seven and a half digits on the way.

7 figures · Block preconditioning rung 4