Exact arithmetic, and what it costs instead

An answer with no error in it

An integer matrix eliminated over the rationals rounds nothing, so the forward error is zero, the residual is the zero vector, and the identity this site is built on has no terms left. The cost does not vanish with the error. It moves into the length of the numbers, where three correct routes differ by four orders of magnitude.

Worth reading first: The condition number is an amplifier · Elimination is a sequence of choices · The exact answer to a nearby problem.

Every field before this one is about an arithmetic that cannot represent its own answers. A float has fifty-three significand bits, the exact result of almost any operation on two of them does not fit, and the whole apparatus of this site — a backward error saying what the algorithm did, a condition number saying what the problem did to it, and their product being all anybody sees — exists to account for the difference.

There is a class of problems where none of that applies. If a matrix has integer entries and the question is its determinant, its rank, or the solution of a system it defines, then every quantity involved is a rational number, every arithmetic operation on rationals produces a rational, and an elimination performed in them rounds nothing at any step. The forward error is zero. The backward error is zero. The residual is the zero vector, not a small vector. The condition number is a perfectly good number and it multiplies nothing.

Three exact eliminations of the same random matrix, and the widest number each one formsEvery route here returns the same determinant exactly, so there is no error to plot and the y-axis is the length of the intermediates in bits. At n = 12 the answer itself is 40 bits and Hadamard's bound allows 52. Fraction-free elimination never forms a number wider than 40 bits, because each of its intermediates is a minor of the original matrix. Rational elimination that reduces every fraction to lowest terms reaches 39 bits and pays a gcd on every arithmetic operation to stay there. Rational elimination that does not reduce reaches 1422165 bits — 35600 times the width of the fraction-free route and 35600 times the width of the answer. All three are correct.3456789101112110¹10²10³10⁴10⁵10⁶nwidest intermediate, in bitsrationals, not reducedfraction-free · rationals reduced · the answerone answer, three widthsn12answer40Hadamard bound52fraction-free40reduced rationals39unreduced1.4·10⁶the error is zero on every curvethe cost is the length of the numbers
Fig. 1 Three exact eliminations of one twelve by twelve integer matrix, plotted against the length of the widest number each one forms. Every route returns the same determinant, so there is no error to draw and the vertical axis is bits.

So the interesting question is not whether it works. It is what it costs, and the answer is the subject of this field: the cost stops being a count of operations and becomes the length of the numbers those operations are performed on.

The identity, with nothing in it

The relation the rest of this site reports every measurement against is

forward error  ⪅  condition number × backward error

and the essay that introduced it is careful that both factors are separately measurable. Here the right-hand side is zero because its second factor is, and a product with a zero in it carries no information about the other one. The matrix above has a two-norm condition number of 32.7. That is a true statement about the matrix and it predicts nothing at all about the computed answer, because the computed answer is the answer.

This is not a rhetorical point. It changes what a figure in this field can be. There is no residual badge to print, because ‖PA − LU‖ is a matrix of exact zeros at every size and on every family, and a badge reading 0 on nine consecutive figures tells a reader nothing they did not know from the first. What can be measured is on the vertical axis of the hero: the number of bits in the widest integer the routine had to hold.

Where the cost went

Take a twelve by twelve matrix whose entries are integers between −9 and 9. Its determinant is 2,269,932,161,262 — twelve digits, forty bits, and small enough to write on a line. Three routes reach it:

  • Fraction-free elimination, which stays in the integers by dividing each step by the previous pivot. The widest number it ever forms is 40 bits.
  • Elimination over the rationals with every fraction reduced to lowest terms, which reaches 39 bits and calls a greatest common divisor on every entry it forms to stay there.
  • Elimination over the rationals without reducing, which is the routine anybody writes first: numerator times numerator, denominator times denominator, subtract with a common denominator. It reaches 1,422,165 bits.
12 × 12: three exact determinants agree to the digit, and the float one does notdet A = 647452350646 on a random matrix of order 12, computed three ways with no rounding anywhere: fraction-free elimination, elimination over the rationals, and elimination modulo 3 word-sized primes put back together by the Chinese remainder theorem. The three differ by exactly zero — as integers, not to a tolerance. The floating-point LU of the same matrix returns 6.47452351·10¹¹, a relative error of 3.394·10⁻¹⁵, which is what a determinant of 40 bits does to a significand of 53. The modular route's entries never exceed 26 bits, which is why the number of primes has to be decided in advance: Hadamard's bound says 53 bits are needed and 3 primes supply them.Bareiss40 bitsrationals39 bitsmodular26 bitsthe answer40 bitsHadamard53 bitswidest number formedrelative error against the exact answerfraction-free (Bareiss)0rational elimination0modular, 3 primes0floating-point LU3.4·10⁻¹⁵det A = 647452350646three routes, one integerand no tolerance anywhere in the comparison
Fig. 2 The same determinant computed four ways. The three exact routes differ by exactly zero, as integers rather than to a tolerance; the floating-point LU differs in the last digits because the answer needs forty bits of significand and it has fifty-three of them holding a product of twelve pivots.

A million and a half bits is four hundred and twenty-eight thousand decimal digits, carried in order to compute a twelve-digit answer. Nothing about that is a rounding error, a stability question or a choice of pivot. All three routes are correct. All three return the identical integer. They differ by a factor of thirty-five thousand in the size of the objects they manipulate, and the flop count — which is the same n³/3 for all three — sees none of it, exactly as the flop count fails to predict the time for a different reason two fields over.

That is the field’s organising fact, and it can be stated in one line: an operation on a forty-thousand-bit integer is not one operation.

Why the unreduced route explodes, exactly

The mechanism is worth writing out because it is short and because it is the reason the other two routes exist.

At each elimination step, an entry becomes aᵢⱼ − (aᵢₖ/aₖₖ)·aₖⱼ. Written over a common denominator with no cancellation performed, the new numerator is a product of two previous numerators minus another such product, and the new denominator is a product of two previous denominators. So if every entry has b bits before the step, every entry has about 2b after it. Twelve steps of doubling from four bits is 4·2¹² = 16,384 — and the measured figure is larger still, because the subtraction of two products of different denominators multiplies them.

Nothing cancels unless something is asked to cancel. The fractions are in lowest terms after the last step, in the sense that the determinant comes out as an integer over one; they are simply never put there on the way.

Three exact eliminations of the same random matrix, and the widest number each one formsEvery route here returns the same determinant exactly, so there is no error to plot and the y-axis is the length of the intermediates in bits. At n = 8 the answer itself is 24 bits and Hadamard's bound allows 32. Fraction-free elimination never forms a number wider than 24 bits, because each of its intermediates is a minor of the original matrix. Rational elimination that reduces every fraction to lowest terms reaches 23 bits and pays a gcd on every arithmetic operation to stay there. Rational elimination that does not reduce reaches 4938 bits — 206 times the width of the fraction-free route and 206 times the width of the answer. All three are correct.345678110¹10²10³10⁴nwidest intermediate, in bitsrationals, not reducedfraction-free · rationals reduced · the answerone answer, three widthsn8answer24Hadamard bound32fraction-free24reduced rationals23unreduced4938the error is zero on every curvethe cost is the length of the numbers
Fig. 3 The same measurement stopped at n = 8, where the unreduced route is at 4,938 bits and the other two are at 24 and 23. The separation is already three orders of magnitude and the curves have not begun to be dramatic.

So the obvious repair is to reduce, and it works: the reduced route’s widest intermediate is 39 bits at n = 12, one bit narrower than the fraction-free route’s. What it costs is that every addition, every multiplication and every subtraction is followed by a greatest common divisor of two integers — and a gcd is not a constant-time operation any more than a multiplication of two long integers is. The reduction is exactly the work the unreduced route was avoiding.

That fork — pay in gcds or pay in length — is the first of the field’s three answers, and the second is the one that makes it unnecessary.

The division that is always exact

Fraction-free elimination performs

aᵢⱼ  ←  (aₖₖ·aᵢⱼ − aᵢₖ·aₖⱼ) / (previous pivot)

and never leaves the integers, because that division is exact at every step of every matrix. It is not exact by luck and it is not exact because of anything about these entries. It is exact because the numerator, before the division, is a determinant with the previous pivot as a factor — which is Sylvester’s identity, and which the next essay in this field checks against the minors themselves rather than quoting.

The consequence for cost is immediate. Each intermediate entry is a minor of the original matrix, so its size is bounded by Hadamard’s bound on a minor of that size, which grows like k·(entry width + log k) — linearly in the step index, not exponentially. That is the flat lower group of curves in the hero, and it is why the fraction-free route needs no gcds: the division that keeps the numbers small is a single exact division rather than a search for a common factor.

Every boxed entry is a 3 × 3 minor of the matrix beside itA 5 × 5 random matrix and its state after 2 steps of fraction-free elimination. Each of the 9 boxed entries was formed as (pivot × entry − column × row) divided by the previous pivot, and every one of those 25 divisions came out exact. The reason is that the entry at (i, j) is exactly the 3 × 3 minor of the ORIGINAL matrix on rows 0…1 and i, columns 0…1 and j — checked here against a cofactor expansion of A rather than against a second elimination. The widest number anywhere in the picture is 9 bits.1-2-6-7-474-26-12-5-8-524-394151435A, the matrix as given1-2-6-7-401840552700112217207003943011710017279153after 2 fraction-free steps3 × 3 minors25 divisions, all exactthe intermediates are minorswhich is why the divisions come out whole
Fig. 4 A five by five matrix and its state after two fraction-free steps. Every boxed entry is a three by three minor of the matrix on the left, checked against a cofactor expansion rather than against a second elimination, and every division that produced it came out whole.

How large the answer is allowed to be

The third route removes growth entirely by removing the room for it: work modulo a prime that fits in a machine word, where no entry can exceed the modulus, and put the answer back together from several such computations by the Chinese remainder theorem.

That only works if the number of primes is known in advance, and it is, because Hadamard’s bound says how large a determinant can be before any of it is computed: |det A| is at most the product of the row two-norms. On the matrix above the bound is 52 bits and the answer occupies 40, so two twenty-six-bit primes suffice and the routine uses two.

The determinant's length, and the bound that decides how many primes to use before any of them is chosenSolid: Hadamard's bound on |det A|, the product of the row 2-norms, in bits. Dotted: the determinant's actual length. The bound is what a modular determinant budgets against — at n = 10 the random family needs 42 bits, which is 2 word-sized primes at 26 bits each, and that count is fixed before the first residue is computed. On a Sylvester Hadamard matrix of order 8 the bound is 13 bits and the determinant has 13, so the bound is attained rather than merely satisfied. On the unit triangular family the determinant is 1 at every size while the bound grows, which is the gap the budget pays for: a matrix can have a determinant of one and rows that are long.345678910110¹10²nbitsthe budget and what it buysrandom, bound42random, actual37primes needed2Hadamard n = 8, bound13Hadamard n = 8, actual13the count is decided by a theorembefore any arithmetic happens
Fig. 5 Hadamard’s bound in bits against the determinant’s actual length, on four families. The gap between a solid line and its own dotted one is the waste in the budget — smallest where the entries are large and used, largest on a family whose determinant is 1 at every size.

The bound is not loose by construction, which matters, because a budget derived from a bound that is never attained would be a budget that is always wrong. On a Sylvester Hadamard matrix of order eight the bound is 25 bits and the determinant occupies 25: the bound is reached. On the Pascal matrix, whose entries reach twenty bits at n = 12 and whose determinant is exactly 1, the bound is 149 bits and the answer is one. Both are correct behaviours of the same inequality, and the second is the field’s standing warning that entry growth during an elimination is not predicted by the size of the answer.

What is the same, and what is not

Three things this site has said elsewhere survive into the exact setting unchanged.

Two routes to every number. The habit is stronger here rather than weaker: two floating-point routes to a quantity agree to a relative 10⁻¹⁵ and the essays price the last digit, but two exact routes agree or they do not. There is no tolerance in the comparison and nothing to negotiate. The determinant above is checked three ways and the differences are the integer zero.

A rank is still a decision, and now for a different reason. The essay that made that argument made it about a threshold: a floating-point matrix has singular values and no rank, and the rank is where somebody puts a line. Remove the arithmetic error and the threshold goes away — and the answer still is not a property of the array of numbers, because it depends on the ring the entries are read in. That is later in this field and it is a genuinely separate failure, closer to a rank that depends on the thread count than to a threshold.

An exact answer to the wrong problem is still wrong. The whole apparatus of conditioning survives, because it was never about arithmetic. It is about the map from data to answer, and if the data was measured rather than given then the exactness buys nothing at all: the field’s last essay measures a case where an exact solve and a float solve of the same physically-obtained system differ by less than either differs from the truth.

8 × 8: three exact determinants agree to the digit, and the float one does notdet A = 3458508270104059228258152 on a wide matrix of order 8, computed three ways with no rounding anywhere: fraction-free elimination, elimination over the rationals, and elimination modulo 4 word-sized primes put back together by the Chinese remainder theorem. The three differ by exactly zero — as integers, not to a tolerance. The floating-point LU of the same matrix returns 3.45850827·10²⁴, a relative error of 3.105·10⁻¹⁶, which is what a determinant of 82 bits does to a significand of 53. The modular route's entries never exceed 26 bits, which is why the number of primes has to be decided in advance: Hadamard's bound says 87 bits are needed and 4 primes supply them.Bareiss82 bitsrationals77 bitsmodular26 bitsthe answer82 bitsHadamard87 bitswidest number formedrelative error against the exact answerfraction-free (Bareiss)0rational elimination0modular, 4 primes0floating-point LU3.1·10⁻¹⁶det A = 3458508270104059228258152three routes, one integerand no tolerance anywhere in the comparison
Fig. 6 The same three routes on a family with three-digit entries, where the determinant occupies 82 bits and the floating-point LU no longer has the significand to hold it. Four primes are budgeted rather than two.

What a float does with a forty-bit answer

The comparison in that figure is worth reading carefully, because it is the only place in this field where floating point appears at all.

A determinant of 2,269,932,161,262 needs 40 bits. A double has 53, so the answer fits — and the float LU still misses it, by a relative 10⁻¹⁵, because the answer is not what is being held. What is being held is a running product of twelve pivots, each of which is a quotient of quantities that do not fit, and the last digits of that product are the ones the cancellation took. Raise the entries to three digits and the determinant needs 82 bits, at which point no amount of care in the elimination can help, because the destination has 53.

That is the honest boundary of the field. Floating point is not bad at this; it is answering a different question, at a fixed precision, for a fixed price. Exact arithmetic answers the asked question at a price that depends on the answer’s own size — and once that size passes the significand, the two are not competing.

The measurement that changed how this is stated

The first version of the swell figure plotted the reduced-rational route as the villain, on the strength of the phrase intermediate expression swell, which is what the literature calls this. The measurement refused it. The reduced route sits at 39 bits at n = 12, one bit below the fraction-free route, and at 13 bits at n = 6 where fraction-free is at 16 — narrower than the answer’s own length, because reducing a fraction whose value is small can throw away a great deal.

The swell belongs to the elimination without the gcd, and the two are separate routines that get called by the same name. Stating it as a property of “rational arithmetic” would have been a claim that a five-line change to the inner loop falsifies. What the measurement supports is the fork: the cost is length or it is gcds, and fraction-free elimination is the third road that pays neither.

Three exact eliminations of the same wide matrix, and the widest number each one formsEvery route here returns the same determinant exactly, so there is no error to plot and the y-axis is the length of the intermediates in bits. At n = 10 the answer itself is 103 bits and Hadamard's bound allows 110. Fraction-free elimination never forms a number wider than 103 bits, because each of its intermediates is a minor of the original matrix. Rational elimination that reduces every fraction to lowest terms reaches 103 bits and pays a gcd on every arithmetic operation to stay there. Rational elimination that does not reduce reaches 754790 bits — 7330 times the width of the fraction-free route and 7330 times the width of the answer. All three are correct.345678910110¹10²10³10⁴10⁵10⁶nwidest intermediate, in bitsrationals, not reducedfraction-free · rationals reduced · the answerone answer, three widthsn10answer103Hadamard bound110fraction-free103reduced rationals103unreduced7.5·10⁵the error is zero on every curvethe cost is the length of the numbers
Fig. 7 On three-digit entries the same picture, shifted up. The lower group starts higher, because the answer is longer, and the separation between the groups is unchanged, because the doubling per step does not know how large the entries were.

What this field is for

Four things are true here that are true nowhere else on this site, and they are why the field exists rather than being a paragraph in the elimination one.

The answer is a certificate. Not a value with an error bar and not a value with a residual — an integer, which either is or is not the determinant, and whose correctness is a statement about arithmetic rather than about a tolerance. The interval-arithmetic essays reach a proved enclosure of a real answer; this reaches the answer.

The cost model is a different one. Everything else on this site counts operations and then argues about where the time actually went — cache, bandwidth, messages. Here the flop count is wrong before any of that, because the operands have lengths, and the length is decided by the algorithm rather than by the machine.

The failure mode is size, not accuracy. A routine here does not return a bad answer; it returns the right one after allocating a gigabyte, or it does not return. That is a different thing to engineer against, and the engineering is Hadamard’s bound and the modular route.

And the ground truth this site has used all along lives here. The Hilbert inverse in closed form, the spanning-tree count from a determinant, the effective resistance as a ratio of two integers: each of those is a float measurement checked against an exact one, and the exact one was computed by the routines this field is about. The instrument has been in use since the collection’s first essays. This is the field that measures the instrument.

At other settings

Three exact eliminations of the same pascal matrix, and the widest number each one formsEvery route here returns the same determinant exactly, so there is no error to plot and the y-axis is the length of the intermediates in bits. At n = 11 the answer itself is 1 bits and Hadamard's bound allows 124. Fraction-free elimination never forms a number wider than 18 bits, because each of its intermediates is a minor of the original matrix. Rational elimination that reduces every fraction to lowest terms reaches 18 bits and pays a gcd on every arithmetic operation to stay there. Rational elimination that does not reduce reaches 18 bits — 1 times the width of the fraction-free route and 18 times the width of the answer. All three are correct.34567891011110¹nwidest intermediate, in bitsrationals, not reducedfraction-free · rationals reduced · the answerone answer, three widthsn11answer1Hadamard bound124fraction-free18reduced rationals18unreduced18the error is zero on every curvethe cost is the length of the numbers
Fig. 8 The Pascal family, whose determinant is 1 at every size while its entries reach twenty bits. The lower group is flat at one bit and Hadamard’s bound climbs away from it.
Three exact eliminations of the same unitriangular matrix, and the widest number each one formsEvery route here returns the same determinant exactly, so there is no error to plot and the y-axis is the length of the intermediates in bits. At n = 12 the answer itself is 1 bits and Hadamard's bound allows 16. Fraction-free elimination never forms a number wider than 1 bits, because each of its intermediates is a minor of the original matrix. Rational elimination that reduces every fraction to lowest terms reaches 1 bits and pays a gcd on every arithmetic operation to stay there. Rational elimination that does not reduce reaches 1 bits — 1 times the width of the fraction-free route and 1 times the width of the answer. All three are correct.34567891011121nwidest intermediate, in bitsrationals, not reducedfraction-free · rationals reduced · the answerone answer, three widthsn12answer1Hadamard bound16fraction-free1reduced rationals1unreduced1the error is zero on every curvethe cost is the length of the numbers
Fig. 9 Unit lower triangular: the elimination has nothing to do, so every curve is flat and the figure is the floor any growth is measured against.
10 × 10: three exact determinants agree to the digit, and the float one does notdet A = 1 on a pascal matrix of order 10, computed three ways with no rounding anywhere: fraction-free elimination, elimination over the rationals, and elimination modulo 4 word-sized primes put back together by the Chinese remainder theorem. The three differ by exactly zero — as integers, not to a tolerance. The floating-point LU of the same matrix returns 0.999999996, a relative error of 4.195·10⁻⁹, which is what a determinant of 1 bits does to a significand of 53. The modular route's entries never exceed 26 bits, which is why the number of primes has to be decided in advance: Hadamard's bound says 101 bits are needed and 4 primes supply them.Bareiss16 bitsrationals16 bitsmodular26 bitsthe answer1 bitsHadamard101 bitswidest number formedrelative error against the exact answerfraction-free (Bareiss)0rational elimination0modular, 4 primes0floating-point LU4.2·10⁻⁹det A = 1three routes, one integerand no tolerance anywhere in the comparison
Fig. 10 A determinant of exactly 1, budgeted at 106 bits because the rows are long.
4 × 4: three exact determinants agree to the digit, and the float one does notdet A = -4386 on a random matrix of order 4, computed three ways with no rounding anywhere: fraction-free elimination, elimination over the rationals, and elimination modulo 1 word-sized prime put back together by the Chinese remainder theorem. The three differ by exactly zero — as integers, not to a tolerance. The floating-point LU of the same matrix returns -4386, a relative error of 0, which is what a determinant of 13 bits does to a significand of 53. The modular route's entries never exceed 26 bits, which is why the number of primes has to be decided in advance: Hadamard's bound says 15 bits are needed and 1 prime supplies them.Bareiss13 bitsrationals12 bitsmodular26 bitsthe answer13 bitsHadamard15 bitswidest number formedrelative error against the exact answerfraction-free (Bareiss)0rational elimination0modular, 1 prime0floating-point LU0det A = -4386three routes, one integerand no tolerance anywhere in the comparison
Fig. 11 Small enough that the float route is exact too, which is the only regime in which the four bars agree.
Every boxed entry is a 2 × 2 minor of the matrix beside itA 4 × 4 random matrix and its state after 1 step of fraction-free elimination. Each of the 9 boxed entries was formed as (pivot × entry − column × row) divided by the previous pivot, and every one of those 9 divisions came out exact. The reason is that the entry at (i, j) is exactly the 2 × 2 minor of the ORIGINAL matrix on rows 0…0 and i, columns 0…0 and j — checked here against a cofactor expansion of A rather than against a second elimination. The widest number anywhere in the picture is 6 bits.1-2-6-7-474-26-12-5-8-524A, the matrix as given1-2-6-70-1-20-3001138370-21-46-52after 1 fraction-free step2 × 2 minors9 divisions, all exactthe intermediates are minorswhich is why the divisions come out whole
Fig. 12 One step, on a four by four, where every boxed entry is an ordinary two by two determinant.
Every boxed entry is a 4 × 4 minor of the matrix beside itA 5 × 5 pascal matrix and its state after 3 steps of fraction-free elimination. Each of the 4 boxed entries was formed as (pivot × entry − column × row) divided by the previous pivot, and every one of those 29 divisions came out exact. The reason is that the entry at (i, j) is exactly the 4 × 4 minor of the ORIGINAL matrix on rows 0…2 and i, columns 0…2 and j — checked here against a cofactor expansion of A rather than against a second elimination. The widest number anywhere in the picture is 5 bits.111111234513610151410203515153570A, the matrix as given11111012340013600014000417after 3 fraction-free steps4 × 4 minors29 divisions, all exactthe intermediates are minorswhich is why the divisions come out whole
Fig. 13 Three steps on Pascal, where the minors are all 1 and the elimination visibly does no work.
The determinant's length, and the bound that decides how many primes to use before any of them is chosenSolid: Hadamard's bound on |det A|, the product of the row 2-norms, in bits. Dotted: the determinant's actual length. The bound is what a modular determinant budgets against — at n = 13 the random family needs 56 bits, which is 3 word-sized primes at 26 bits each, and that count is fixed before the first residue is computed. On a Sylvester Hadamard matrix of order 8 the bound is 13 bits and the determinant has 13, so the bound is attained rather than merely satisfied. On the unit triangular family the determinant is 1 at every size while the bound grows, which is the gap the budget pays for: a matrix can have a determinant of one and rows that are long.345678910111213110¹10²nbitsthe budget and what it buysrandom, bound56random, actual44primes needed3Hadamard n = 8, bound13Hadamard n = 8, actual13the count is decided by a theorembefore any arithmetic happens
Fig. 14 The budget out to n = 13, where the wide family needs six primes and the random one three.
The exact solution of a 13×13 Hilbert system beside the computed oneTwo columns of numbers: the exact answer, which is the integers one to thirteen, and the answer double-precision elimination returns, with the number of correct digits beside each.H13 x = b, b formed exactly so that x = (1, 2, …, 13)12345678910111213exact1.00002.00003.00133.97865.18864.993010.46720.048921.2657-2.575319.21478.906013.5113computed6.6 correct digits4.8 correct digits3.4 correct digits2.3 correct digits1.4 correct digit0.8 correct digitno correct digitsno correct digitsno correct digitsno correct digitsno correct digits0.6 correct digit1.4 correct digitbackward error2.2·10⁻¹⁷κ = 1.7·10¹⁸The algorithm solved a neighbouring problem perfectly. That problem's answer is this one.right-hand side built in BigInt rationalsthe truth is known
Fig. 15 The instrument this field measures: a float Hilbert solve against the rational answer.
Backward error, forward error and the condition number, with measured valuesTwo boxes at the top — the problem posed and the nearby problem the algorithm answered exactly — and two answers below them, with the distances between all four labelled by numbers from a Hilbert solve.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
Fig. 16 The identity with both of its factors, in the field where they are not zero.
Backward and forward error against the condition numberA log–log plot over twelve decades of condition number. The backward error is a flat line at ten to the minus sixteen; the forward error rises in proportion to the condition number.110²10⁴10⁶10⁸10¹⁰10¹²10¹⁴10⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹condition number κ(A)relative errorforward errorbackward errorpredicted: κ · u8×8, 20 seeds per κ; dashed is the worstthe problem worsens, not the method
Fig. 17 A small residual and a large forward error, which is the situation this field cannot be in.
How much a perturbation of the right-hand side is amplified, κ = 10⁶The cumulative distribution of the amplification factor over two hundred random perturbation directions, with the condition number marked as the upper limit.110¹10²10³10⁴10⁵10⁶10⁷00.250.50.751amplification of the input perturbationfraction of directions at or belowκ = 10·10⁵worst found 7.6·10⁵6×6, 200 directionsmedian reaches 0.29 of κ
Fig. 18 The amplifier that belongs to the problem, and survives the arithmetic being exact.
Growth factor under partial pivoting to n = 40: the bound, the worst case, and realityGrowth factor against matrix size on a logarithmic vertical axis. The two-to-the-n bound rises as a straight line; Wilkinson's matrix sits exactly on it; random matrices stay near one.0816243240110²10⁴10⁶10⁸10¹⁰10¹²10¹⁴matrix size ngrowth factor max|u| / max|a|the 2ⁿ⁻¹ boundworst of 30 randommedian randomWilkinson's matrix sits on the bound30 Gaussian matrices per sizeat n = 40: bound 5.5·10¹¹, worst 4.8
Fig. 19 The elimination field’s own bound that is never attained, against Hadamard’s, which is.
Forward error of a 6×6 Hilbert solve at eight precisionsA bar for each significand width from 12 to 53 bits showing the relative error in the computed solution, with the condition number times the unit roundoff marked as a prediction.κ = 1.5·10⁷ · the exact answer is (1, 2, …, 6)12 bits3.316 bits0.8620 bits4.324 bits0.09330 bits7.6·10⁻⁴36 bits6.2·10⁻⁶43 bits2·10⁻⁷53 bits4.5·10⁻¹¹dashed: κ · unit roundoffone matrix, eight arithmeticsmeasured against a known answer
Fig. 20 The site’s precision knob, which this field has removed.
Relative error of two algebraically identical expressions for (1 − cos x)/x², in binary32A log–log plot of relative error against x. The expression written as it reads loses accuracy as x falls and is eventually wrong in every digit; the rearranged form stays at rounding level.10⁻⁸10⁻⁶10⁻⁴10⁻²110⁻⁹10⁻⁶10⁻³1xrelative error of the computed value(1 − cos x)/x², as written2 sin²(x/2)/x²no digits left at allbinary32 throughoutone function, two spellings · zero below 3.5·10⁻⁴
Fig. 21 The cancellation that takes the last digits of a float determinant, in the field that named it.
Singular values of a rank-4 matrix with noise of relative size 10⁻⁸Ten singular values on a logarithmic axis. The first four sit near one; the rest sit at the noise level, and the vertical distance between the two groups is the evidence for the rank.1234567891010⁻¹⁸10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1indexsingular valuecutoff, σ₁ · 10⁻¹⁰numerical rank 10gap 8.2·10⁶an opiniontrue rank 410×10, built with 4 nonzero valuesrank is a decision
Fig. 22 Rank as a threshold, which is the version of the question this field replaces.
The 12×12 grid Laplacian and its Cholesky factor, ordered by naturalTwo square sparsity plots side by side. The left shows the nonzeros of the matrix; the right shows the nonzeros of its Cholesky factor, with the entries created by elimination marked in a second colour.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
Fig. 23 Fill in a sparse elimination, which is the other way an intermediate becomes larger than its input.
Three exact eliminations of the same random matrix, and the widest number each one formsEvery route here returns the same determinant exactly, so there is no error to plot and the y-axis is the length of the intermediates in bits. At n = 7 the answer itself is 24 bits and Hadamard's bound allows 27. Fraction-free elimination never forms a number wider than 24 bits, because each of its intermediates is a minor of the original matrix. Rational elimination that reduces every fraction to lowest terms reaches 17 bits and pays a gcd on every arithmetic operation to stay there. Rational elimination that does not reduce reaches 711 bits — 29.6 times the width of the fraction-free route and 29.6 times the width of the answer. All three are correct.34567110¹10²10³nwidest intermediate, in bitsrationals, not reducedfraction-free · rationals reduced · the answerone answer, three widthsn7answer24Hadamard bound27fraction-free24reduced rationals17unreduced711the error is zero on every curvethe cost is the length of the numbers
Fig. 24 Seven by seven, where the unreduced route is already at 711 bits and the others at 24.
Three exact eliminations of the same wide matrix, and the widest number each one formsEvery route here returns the same determinant exactly, so there is no error to plot and the y-axis is the length of the intermediates in bits. At n = 12 the answer itself is 121 bits and Hadamard's bound allows 133. Fraction-free elimination never forms a number wider than 121 bits, because each of its intermediates is a minor of the original matrix. Rational elimination that reduces every fraction to lowest terms reaches 119 bits and pays a gcd on every arithmetic operation to stay there. Rational elimination that does not reduce reaches 2926485 bits — 24200 times the width of the fraction-free route and 24200 times the width of the answer. All three are correct.3456789101112110¹10²10³10⁴10⁵10⁶nwidest intermediate, in bitsrationals, not reducedfraction-free · rationals reduced · the answerone answer, three widthsn12answer121Hadamard bound133fraction-free121reduced rationals119unreduced2.9·10⁶the error is zero on every curvethe cost is the length of the numbers
Fig. 25 Three-digit entries out to twelve, where the answer needs 121 bits and the budget six primes.
8 × 8: three exact determinants agree to the digit, and the float one does notdet A = 1 on a unitriangular matrix of order 8, computed three ways with no rounding anywhere: fraction-free elimination, elimination over the rationals, and elimination modulo 1 word-sized prime put back together by the Chinese remainder theorem. The three differ by exactly zero — as integers, not to a tolerance. The floating-point LU of the same matrix returns 1, a relative error of 0, which is what a determinant of 1 bits does to a significand of 53. The modular route's entries never exceed 26 bits, which is why the number of primes has to be decided in advance: Hadamard's bound says 10 bits are needed and 1 prime supplies them.Bareiss1 bitsrationals1 bitsmodular26 bitsthe answer1 bitsHadamard10 bitswidest number formedrelative error against the exact answerfraction-free (Bareiss)0rational elimination0modular, 1 prime0floating-point LU0det A = 1three routes, one integerand no tolerance anywhere in the comparison
Fig. 26 A determinant of exactly 1, reached by every route including the float one.
6 × 6: three exact determinants agree to the digit, and the float one does notdet A = 2435091120 on a hilbert matrix of order 6, computed three ways with no rounding anywhere: fraction-free elimination, elimination over the rationals, and elimination modulo 4 word-sized primes put back together by the Chinese remainder theorem. The three differ by exactly zero — as integers, not to a tolerance. The floating-point LU of the same matrix returns 2.43509112·10⁹, a relative error of 1.329·10⁻¹¹, which is what a determinant of 32 bits does to a significand of 53. The modular route's entries never exceed 26 bits, which is why the number of primes has to be decided in advance: Hadamard's bound says 85 bits are needed and 4 primes supply them.Bareiss40 bitsrationals15 bitsmodular26 bitsthe answer32 bitsHadamard85 bitswidest number formedrelative error against the exact answerfraction-free (Bareiss)0rational elimination0modular, 4 primes0floating-point LU1.3·10⁻¹¹det A = 2435091120three routes, one integerand no tolerance anywhere in the comparison
Fig. 27 The integer Hilbert matrix, where conditioning costs the float route and nothing else.
The determinant's length, and the bound that decides how many primes to use before any of them is chosenSolid: Hadamard's bound on |det A|, the product of the row 2-norms, in bits. Dotted: the determinant's actual length. The bound is what a modular determinant budgets against — at n = 7 the random family needs 27 bits, which is 2 word-sized primes at 26 bits each, and that count is fixed before the first residue is computed. On a Sylvester Hadamard matrix of order 8 the bound is 13 bits and the determinant has 13, so the bound is attained rather than merely satisfied. On the unit triangular family the determinant is 1 at every size while the bound grows, which is the gap the budget pays for: a matrix can have a determinant of one and rows that are long.34567110¹10²nbitsthe budget and what it buysrandom, bound27random, actual24primes needed2Hadamard n = 8, bound13Hadamard n = 8, actual13the count is decided by a theorembefore any arithmetic happens
Fig. 28 The budget at small sizes, where one prime does for every family.
Every boxed entry is a 3 × 3 minor of the matrix beside itA 6 × 6 random matrix and its state after 2 steps of fraction-free elimination. Each of the 16 boxed entries was formed as (pivot × entry − column × row) divided by the previous pivot, and every one of those 41 divisions came out exact. The reason is that the entry at (i, j) is exactly the 3 × 3 minor of the ORIGINAL matrix on rows 0…1 and i, columns 0…1 and j — checked here against a cofactor expansion of A rather than against a second elimination. The widest number anywhere in the picture is 9 bits.1-2-6-7-474-26-12-5-8-524-394151435-6-5-33-29-46184A, the matrix as given1-2-6-7-4706302718-3300354255168-30300-96-69-4214700308466-9000-60612108after 2 fraction-free steps3 × 3 minors41 divisions, all exactthe intermediates are minorswhich is why the divisions come out whole
Fig. 29 Six by six after two steps, where every boxed entry is a three by three minor.
The exact solution of a random system is longer than the systemAx = b with A a random integer matrix and b chosen so the answer is a genuine vector of fractions. At n = 9 the question occupies 324 bits — 9² entries of 4 bits each — and the answer occupies 454, a ratio of 1.4. Every entry is a ratio of two determinants by Cramer's rule, checked here by multiplying each one back by det A and comparing it with the determinant of the column-replaced matrix as an integer — so every denominator divides det A and none of them is longer than its 25 bits, while a numerator can be, at 28 here. The answer IS that long, whatever route produced it. That is a floor no exact solver can get under, and it is why the modular route reconstructs the answer rather than eliminating its way to it.345678910²nbitsthe answerthe question · det Abits in, bits outn9question324answer454det A25widest part28the answer is the floorand no route writes it down more cheaply
Fig. 30 The other floor: how long the answer to a system of this kind is.
5 × 5: the invariant factors, computed twice by routes with no algorithm in commonThe Smith normal form of a random matrix of order 5 is the diagonal [1, 1, 1, 1, 22616], reached by unimodular row and column operations, with each entry dividing the next — a property the elimination does not give and which has to be enforced and asserted. The same numbers computed from the definition — the gcd of all k × k minors, over the gcd of all (k−1) × (k−1) minors — are [1, 1, 1, 1, 22616], and the two agree entry for entry. Their product is 22616, which is |det A| = 22616. The Hermite form of the same matrix reaches 24 bits against the Smith form's 15, and its transform has determinant -1, so the lattice the rows generate is unchanged.s10 bitss20 bitss30 bitss40 bitss514 bitsinvariant factors, in bitstwo routes, and a conserved quantitydet A-2.3·10⁴Π invariants2.3·10⁴SNF widest15HNF widest24det of the transform-1an algorithm and a definitionagreeing as integers
Fig. 31 The invariant factors of a matrix of this family, which its determinant is the product of.
det A = 30030, and its residue at 36 primes — 5 of which divide itEach dot is det A mod p, drawn as a fraction of p so the primes are comparable. A dot on the floor is a prime at which the elimination reports a singular matrix — and the matrix is not singular: its determinant is 30030. 5 of 36 primes here are unlucky, a rate of 13.9 per cent, which is what it looks like when the primes are small. At the 26-bit primes a modular determinant actually uses, a prime divides a 15-bit determinant with probability of order 2.24·10⁻⁷ — but the consequence is unchanged, and there is no way to see it from one residue. What sees it is the budget: computing more residues than Hadamard's bound requires, and refusing to reconstruct when the extra ones disagree.det mod p, as a fraction of p311a prime that divides the answerdet A3·10⁴primes swept36unlucky5rate0.14det, in bits15singular mod p is not singularand one residue cannot tell them apart
Fig. 32 The risk the modular route runs in exchange for having no growth.
An exact solve of a measured 8 × 8 system, against the error the measurement put thereThe rising curve is the distance from the exact solution of a perturbed system to the exact solution of the unperturbed one, as the perturbation grows from a relative 10⁻¹⁶ to 1. Every point on it was computed with no rounding anywhere: the elimination is over the rationals, the residual of each solve is the zero vector, and the error is entirely the data's. At a measurement accurate to 10^-8 the answer is wrong by 2.049·10⁻⁸, which is 2.05 times the data's own error — the matrix's condition number, doing the only job left for it to do. The straight line of slope one is the whole finding: removing the arithmetic error changes nothing about this quantity, because this quantity was never about the arithmetic.10⁻¹⁶10⁻¹⁴10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²110⁻¹⁶10⁻¹⁴10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²1relative error in the datarelative error in the answerthe answer's errorthe data's errornothing here is the arithmetic'sdata error10⁻⁸answer error2·10⁻⁸amplification2rounding committed0residual0the residual is exactly zeroand the answer is wrong anyway
Fig. 33 A well-conditioned system, where the exactness and the data agree about the answer.
Three exact eliminations of the same random matrix, and the widest number each one formsEvery route here returns the same determinant exactly, so there is no error to plot and the y-axis is the length of the intermediates in bits. At n = 5 the answer itself is 15 bits and Hadamard's bound allows 18. Fraction-free elimination never forms a number wider than 15 bits, because each of its intermediates is a minor of the original matrix. Rational elimination that reduces every fraction to lowest terms reaches 15 bits and pays a gcd on every arithmetic operation to stay there. Rational elimination that does not reduce reaches 85 bits — 5.67 times the width of the fraction-free route and 5.67 times the width of the answer. All three are correct.345110¹10²nwidest intermediate, in bitsrationals, not reducedfraction-free · rationals reduced · the answerone answer, three widthsn5answer15Hadamard bound18fraction-free15reduced rationals15unreduced85the error is zero on every curvethe cost is the length of the numbers
Fig. 34 Five by five, where the unreduced route first leaves the group.
Three exact eliminations of the same pascal matrix, and the widest number each one formsEvery route here returns the same determinant exactly, so there is no error to plot and the y-axis is the length of the intermediates in bits. At n = 9 the answer itself is 1 bits and Hadamard's bound allows 80. Fraction-free elimination never forms a number wider than 14 bits, because each of its intermediates is a minor of the original matrix. Rational elimination that reduces every fraction to lowest terms reaches 14 bits and pays a gcd on every arithmetic operation to stay there. Rational elimination that does not reduce reaches 14 bits — 1 times the width of the fraction-free route and 14 times the width of the answer. All three are correct.3456789110¹nwidest intermediate, in bitsrationals, not reducedfraction-free · rationals reduced · the answerone answer, three widthsn9answer1Hadamard bound80fraction-free14reduced rationals14unreduced14the error is zero on every curvethe cost is the length of the numbers
Fig. 35 Pascal at nine, a determinant of one under a bound of eighty bits.
6 × 6: three exact determinants agree to the digit, and the float one does notdet A = -522130948717689337 on a wide matrix of order 6, computed three ways with no rounding anywhere: fraction-free elimination, elimination over the rationals, and elimination modulo 3 word-sized primes put back together by the Chinese remainder theorem. The three differ by exactly zero — as integers, not to a tolerance. The floating-point LU of the same matrix returns -5.22130949·10¹⁷, a relative error of 2.451·10⁻¹⁶, which is what a determinant of 59 bits does to a significand of 53. The modular route's entries never exceed 26 bits, which is why the number of primes has to be decided in advance: Hadamard's bound says 65 bits are needed and 3 primes supply them.Bareiss59 bitsrationals59 bitsmodular26 bitsthe answer59 bitsHadamard65 bitswidest number formedrelative error against the exact answerfraction-free (Bareiss)0rational elimination0modular, 3 primes0floating-point LU2.5·10⁻¹⁶det A = -522130948717689337three routes, one integerand no tolerance anywhere in the comparison
Fig. 36 Three-digit entries at order six, budgeted at three primes.
8 × 8: three exact determinants agree to the digit, and the float one does notdet A = 4096 on a hadamard matrix of order 8, computed three ways with no rounding anywhere: fraction-free elimination, elimination over the rationals, and elimination modulo 1 word-sized prime put back together by the Chinese remainder theorem. The three differ by exactly zero — as integers, not to a tolerance. The floating-point LU of the same matrix returns 4096, a relative error of 0, which is what a determinant of 13 bits does to a significand of 53. The modular route's entries never exceed 26 bits, which is why the number of primes has to be decided in advance: Hadamard's bound says 14 bits are needed and 1 prime supplies them.Bareiss13 bitsrationals4 bitsmodular26 bitsthe answer13 bitsHadamard14 bitswidest number formedrelative error against the exact answerfraction-free (Bareiss)0rational elimination0modular, 1 prime0floating-point LU0det A = 4096three routes, one integerand no tolerance anywhere in the comparison
Fig. 37 The case where Hadamard’s bound is attained rather than satisfied.
Every boxed entry is a 4 × 4 minor of the matrix beside itA 6 × 6 wide matrix and its state after 3 steps of fraction-free elimination. Each of the 9 boxed entries was formed as (pivot × entry − column × row) divided by the previous pivot, and every one of those 50 divisions came out exact. The reason is that the entry at (i, j) is exactly the 4 × 4 minor of the ORIGINAL matrix on rows 0…2 and i, columns 0…2 and j — checked here against a cofactor expansion of A rather than against a second elimination. The widest number anywhere in the picture is 39 bits.-257325-86896-994-523-996278116-612-894131-704-870-29973890972-390570-7246836-574625-308989514123257-105-573652726-105A, the matrix as given-257325-86896-994-523016241-28664953925887817798100-1.756·10⁸-1.159·10⁶3.545·10⁷1.062·10⁸0009.927·10¹⁰-8.591·10¹⁰-2.483·10¹⁰000-1.421·10¹¹-5.391·10¹¹-4.05·10¹¹000-1.157·10¹¹1.006·10¹¹2.743·10¹¹after 3 fraction-free steps4 × 4 minors50 divisions, all exactthe intermediates are minorswhich is why the divisions come out whole
Fig. 38 Wide entries, three steps in, every division still exact.
The determinant's length, and the bound that decides how many primes to use before any of them is chosenSolid: Hadamard's bound on |det A|, the product of the row 2-norms, in bits. Dotted: the determinant's actual length. The bound is what a modular determinant budgets against — at n = 12 the random family needs 52 bits, which is 3 word-sized primes at 26 bits each, and that count is fixed before the first residue is computed. On a Sylvester Hadamard matrix of order 8 the bound is 13 bits and the determinant has 13, so the bound is attained rather than merely satisfied. On the unit triangular family the determinant is 1 at every size while the bound grows, which is the gap the budget pays for: a matrix can have a determinant of one and rows that are long.3456789101112110¹10²nbitsthe budget and what it buysrandom, bound52random, actual40primes needed3Hadamard n = 8, bound13Hadamard n = 8, actual13the count is decided by a theorembefore any arithmetic happens
Fig. 39 The budget out to twelve, where the families have separated completely.
5 × 5: the invariant factors, computed twice by routes with no algorithm in commonThe Smith normal form of a wide matrix of order 5 is the diagonal [1, 1, 1, 1, 2580458337283640], reached by unimodular row and column operations, with each entry dividing the next — a property the elimination does not give and which has to be enforced and asserted. The same numbers computed from the definition — the gcd of all k × k minors, over the gcd of all (k−1) × (k−1) minors — are [1, 1, 1, 1, 2580458337283640], and the two agree entry for entry. Their product is 2580458337283640, which is |det A| = 2580458337283640. The Hermite form of the same matrix reaches 106 bits against the Smith form's 52, and its transform has determinant -1, so the lattice the rows generate is unchanged.s10 bitss20 bitss30 bitss40 bitss551 bitsinvariant factors, in bitstwo routes, and a conserved quantitydet A-2.6·10¹⁵Π invariants2.6·10¹⁵SNF widest52HNF widest106det of the transform-1an algorithm and a definitionagreeing as integers
Fig. 40 The invariant factors of a wide matrix, one of them fourteen digits.
Recovering 9973/113 from its residue: nothing below 2^28, everything aboveThe fraction 9973/113 is reduced modulo a prime just above 2^e and then reconstructed from that one residue by the extended Euclidean algorithm. Below e = 28 the reconstruction returns a different fraction or refuses; at and above it, it returns 9973/113 exactly, at every e tried. The crossover is not a property of the algorithm: a residue determines a rational uniquely once the modulus exceeds twice the square of the LARGER of the two parts, which here is 2^27.57, and below that there are two fractions with the same residue and both parts small enough. The measured first success is 2^28. There is no band in which the answer is is nearly right — the value returned below the cliff is a perfectly good fraction that is not this one.2^28modulus, as a power of twoa lattice with one short vectornumerator9973denominator1132·max(n, d)², bits28first recovered at28moduli tried42below the bound there are two answersand the algorithm cannot prefer one
Fig. 41 The bound the modular route’s second stage works to.

What links here

Computed from the collection, not written here: the essays that point at this one.

Reads more easily once this is understood

Essays that name this one as worth reading first.

Shares its objects with

Essays that name at least two of the same things, and that neither author linked.

Named objects

A flat tag is an object no other essay names yet.

Backward errorBit lengthCondition numberExact arithmeticExact ground truthFraction-free eliminationHadamard boundIntermediate expression swell