Numerical Linear Algebra
where the algebra stops being the arithmetic
The arithmetic underneath

What a float can hold

The representable numbers are not a fine fuzz spread evenly over the line. They are evenly spaced inside each power-of-two interval and twice as far apart in the next one up, and almost everything else in this subject is a consequence of that one fact.
15 min read 6 figures Where the algebra stops being the arithmeticPrecision as the knobTwo routes to a number

Ask what a computer does to a number and the usual answer is that it stores an approximation. That is true and it is nearly useless, because it leaves the impression of a uniform haze — as if every real number were nudged to the nearest point of a fine, evenly spaced grid. The grid is not evenly spaced. It is a sequence of grids, each twice as coarse as the one below it, joined end to end at the powers of two.

The representable numbers with a 3-bit significandA number line from 0.5 to 4 with a tick at every representable value. The ticks are evenly spaced inside each power-of-two interval and twice as far apart in the next one up.[½, 1)[1, 2)[2, 4)0.5124gap 0.125gap 0.25 — twice as wide8 values per octavespacing doubles at each power of two
Fig. 1 Every value representable with a small significand, between a half and four. Inside each power-of-two interval the values are evenly spaced; crossing into the next interval, the spacing doubles. Drag the slider to add bits: each one halves the gap and doubles the count, and at the fifty-three bits of a double there are 2⁵² values in every octave, which is why a real number line looks solid.

That is the whole structure. A floating-point number is a significand — a fixed number of binary digits — multiplied by a power of two. The significand has a fixed number of digits, not a fixed size of digit, so the absolute precision available scales with the magnitude of the number being stored. Near 1 a double resolves differences of about 2.2·10⁻¹⁶. Near a million it resolves differences of about 1.2·10⁻¹⁰. Near 2⁵³ it resolves differences of exactly 1, and above that it cannot count.

Relative, not absolute

The useful way to hold this is that floating point makes a promise about relative error and no promise at all about absolute error. Store any real number in the normal range and what comes back differs from it by at most a fixed fraction of itself. That fraction is the unit roundoff, written u, and for a double it is 2⁻⁵³ ≈ 1.11·10⁻¹⁶.

Two conventions circulate and they differ by a factor of two, which is a small nuisance worth naming once. Unit roundoff is the largest relative error a single correctly rounded operation can make: half the gap, because rounding goes to the nearer neighbour. Machine epsilon is often defined instead as the gap itself just above 1 — 2⁻⁵² ≈ 2.22·10⁻¹⁶ — which is the smallest number that can be added to 1 and get something other than 1. Both appear in the literature under both names. Everything on this site quotes the gap where a gap is meant and u where a bound is meant, and says which.

The spacing between consecutive numbers at 53-bit precisionA log–log staircase of the gap between neighbouring representable numbers against magnitude. The gap doubles at every power of two and reaches one whole unit partway along.110³10⁶10⁹10¹²10¹⁵10¹⁸10⁻¹⁸10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³110³magnitude of the numbergap to the next representable numbera gap of one whole unitat 1: 2.2·10⁻¹⁶at a million: 1.2·10⁻¹⁰gap reaches 1: 153-bit significandthe gap follows the magnitude
Fig. 2 The gap to the next representable number, plotted against magnitude. It is a staircase, doubling at every power of two, and it crosses one whole unit partway along — which is the reason a loop counting upward by ones through 10¹⁷ never terminates. The generator behind it takes the significand width as a parameter and the whole staircase slides down when it changes, keeping its shape exactly — the slider is turned off here because the page already spends its interaction budget on the number line above, and one draggable figure per page is a rule this site enforces rather than a preference.

The staircase crossing the line at height 1 is worth pausing on. Beyond 2⁵³ the spacing of doubles exceeds one, so consecutive integers are no longer all representable and x + 1 == x becomes true. This is not an obscure boundary: 2⁵³ is about 9·10¹⁵, which a nanosecond timestamp reaches in a little over three months, and which an identifier space reaches routinely. The failure is silent in the worst way — the addition succeeds, returns a number of the right magnitude, and the loop runs forever.

What “the spacing doubles” costs

Every operation on this site inherits the same structure. When two numbers are added, the exact sum is computed and then rounded to the nearest representable value — so the error introduced is at most half the gap at the size of the result. Add two numbers near 1 and the error is around 10⁻¹⁶. Add two numbers near 10⁸ and the error is around 10⁻⁸. Nothing has gone wrong in the second case; the arithmetic is doing exactly what it promised, which was a relative guarantee.

The consequence is that a sum of numbers with very different magnitudes loses the small ones. Take one hundred million and add one to it, a hundred thousand times, in single precision. At 10⁸ the gap between single-precision numbers is 8, so adding 1 rounds straight back to where it started. The answer is 10⁸, the true answer is 100,100,000, and the relative error is 9.99·10⁻⁴. Every one of the hundred thousand additions succeeded.

The order they are added in is about what can be done with that, and the answer is a good deal more than most people expect. But the failure begins here, in the spacing, and it is not repaired by a wider float — a wider float moves the staircase down without changing its shape.

Simulating a precision, and checking the simulation

Everything on this site that says float32 or twenty-four bits is a simulation. The arithmetic is performed in double and each result is rounded to the requested number of significand bits, with ties going to even, which is what IEEE 754 does. That is what makes the precision a parameter that can be dragged rather than a constant baked into the language.

A simulation of arithmetic is exactly the sort of thing that can be subtly wrong in a way no figure would reveal — the pictures would still be smooth, still be monotone, still tell a plausible story. So it is checked against the only independent implementation available: Math.fround, which is the processor’s own single-precision rounding reached from JavaScript. Over twenty thousand test values spread across twenty-four orders of magnitude, the two agree on every single one, bit for bit. That is two routes to a number applied to the ruler rather than to what is being measured, and it is worth doing first.

Two limits are worth stating rather than glossing. The simulation models the significand and not the exponent range, so there is no overflow to infinity at 10³⁸ and no gradual underflow. Those are real behaviours of a real float32 and nothing here claims to show them. And for division and square root the double result is itself rounded before being rounded again, so a doubly rounded value can differ from a singly rounded one in the last bit — for addition, subtraction and multiplication of values at 24 bits or fewer the exact result fits in a double and the question does not arise.

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. 3 One computation at eight precisions. The same Hilbert system, solved by the same code, with the significand width the only thing that changes; the answer is compared against the exact one, and the dashed marks are the prediction κ·u. The error tracks the prediction over eleven orders of magnitude, which is the sense in which precision buys accuracy — a constant factor, bought a bit at a time.

Three numbers worth carrying around

Most of the reasoning in this subject can be done in the head with three quantities, and it is worth having them rather than looking them up.

binary32 binary64
significand bits 24 53
gap just above 1 1.19·10⁻⁷ 2.22·10⁻¹⁶
decimal digits carried about 7 about 16
integers exact up to 1.7·10⁷ 9.0·10¹⁵
√u, which several thresholds sit at 2.4·10⁻⁴ 1.5·10⁻⁸

The last row is the one that surprises people, and it appears three times on this site. A great many numerical thresholds sit not at u but at its square root, because the quantity that fails is a square: ε² falling below the gap next to 1, a condition number squared exceeding what the arithmetic carries, a Newton step whose quadratic term has been rounded away. In single precision √u is 2.4·10⁻⁴, which is not a small number at all — it is a tenth of a percent, and quantities of that size occur in ordinary data.

The road that squares the problem is the sharpest instance: solving least squares through the normal equations fails, in single precision, at exactly that ε, and the failure is not a gentle loss of digits but an exactly singular matrix.

What rounding does to a subtraction

The most-repeated advice about floating point is not to subtract two nearly equal numbers, and the usual explanation for it is wrong in an instructive way. The subtraction is not where the error comes from. If two floating-point numbers are within a factor of two of each other, their difference is exactly representable and the subtraction introduces no error at all — a result due to Sterbenz and easy to see from the staircase, since the difference is a multiple of the common gap.

What the subtraction does is reveal error that was already there. Each operand carries an error of up to half a gap at its own magnitude. When the operands are close and the difference is small, that inherited error is unchanged in absolute terms and enormous relative to the answer.

Relative error of two algebraically identical expressions for (1 − cos x)/x²A log–log plot of relative error against x. The expression written as it reads loses accuracy below x = 10⁻⁴ and is entirely wrong by 10⁻⁸; the rearranged form stays at rounding level.10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²110⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹xrelative error of the computed value(1 − cos x)/x², as written2 sin²(x/2)/x²no digits left at alldouble precision throughoutone function, two spellings
Fig. 4 Two ways of writing the same function, and their error against a series that shares no arithmetic with either. Written as it reads, (1 − cos x)/x² has lost four digits by x = 10⁻⁶ and returns exactly zero from 10⁻⁸ down — not approximately zero, exactly, because 1 − cos x has rounded to nothing and zero divided by anything is zero. The rearranged form holds fifteen digits across the whole range.

Cancellation takes the answer works through what to do about it. The short version is that the fix is never more precision — the naive form fails at 10⁻⁸ in double and would fail at 10⁻¹⁷ in a quadruple — it is a different expression for the same function.

Where the gap becomes a plan

Once the staircase is in view, three things that are usually taught as separate cautions collapse into one.

Do not compare floating-point numbers for equality is not a rule about floating point. It is a rule about the gap: two computations of the same quantity differ by a few units in the last place, and the last place is a different size depending on the magnitude. The right comparison is relative, scaled to the numbers involved, and where the numbers may be near zero it has to be absolute instead — which is why the tolerance helpers used throughout this site take both cases and the tolerance itself is bracketed between a measured noise floor and the smallest deviation that would be a real failure.

Do not subtract nearly equal numbers is the subject of the next essay, and it is the sharpest consequence of the spacing. The subtraction itself is exact — that is the surprise — and the damage was done before it, when the two operands were rounded.

Do not accumulate in a loop without thinking is what the spacing does over many operations, and the interesting part is that thinking helps enormously: the same additions, performed in a different order, are a thousand times more accurate.

What this buys the rest of the site

The reason to start here rather than with matrices is that every matrix claim later on is a claim about accumulated rounding, and accumulated rounding is unreadable without a feel for the size of a single step.

When Gaussian elimination is run without pivoting and returns an answer wrong in its second component, the reason is that a multiplier of 10¹⁷ dragged a useful number up to a magnitude where the thing being added to it had fallen below the gap. When classical Gram–Schmidt returns columns that are not perpendicular, the reason is that a projection coefficient was computed against a vector still containing the component being removed, and the difference between the two coefficients was below the gap and therefore invisible.

Neither of those is a mystery once the staircase is in view. Both are inexplicable without it.

Everything here is reproducible, and that is a decision

A figure about rounding that changed slightly on every build would be worthless, and it would be easy to produce by accident: one call to a system random number generator anywhere in the pipeline and every page acquires a different set of last digits each time it is rendered.

So nothing here calls the system generator. The random matrices are drawn from a counter-based generator written out in thirty lines, with an explicit seed, checked for flatness in two dimensions because that is where a bad generator shows and half these figures plot one computed quantity against another. A figure is byte-identical on every build, and anyone with the seed gets the same matrix.

Where a claim is about typical behaviour rather than about one matrix, it is run across many seeds and the figure reports what held for all of them. One matrix is an anecdote, and a plot of one matrix’s error against its condition number is a plot of a single point pretending to be a law. Both the growth factor and the loss of orthogonality are reported that way here: many seeds, the median drawn solid and the worst of them drawn dashed above it.

Relative error of three summation algorithms in binary32A log–log plot of relative error against the number of terms for naive, pairwise and compensated summation, each measured against the exactly rounded sum.10¹10²10³10⁴10⁵10⁶10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹number of terms addedrelative error against the exact sumin orderin a treecompensatedbinary32 · terms are 1/icompensated: 3·10⁻⁸
Fig. 5 Three summation algorithms in single precision, each measured against the exactly rounded sum of the same terms. Adding a million reciprocals in order loses 2.5·10⁻³ of the answer; adding them in a tree loses 1.9·10⁻⁷; carrying the lost part forward loses 1.2·10⁻⁸. Same numbers, same arithmetic, same count of additions.

The two questions worth carrying

For any computed quantity there are two separate questions, and most confusion in this subject comes from answering one and believing the other has been answered.

The first is: how much error did this arithmetic introduce? That is a question about the algorithm and the precision, and the answer is usually “a few units in the last place per operation, accumulated over however many operations there were”. It is the question this essay is about.

The second is: how much does that error matter to the answer? That is a question about the problem, not about the arithmetic at all, and the answer can be anything from “not at all” to “the answer has no correct digits”. The exact answer to a nearby problem is about separating those two, which is the single most useful idea in the subject and the one this whole site is arranged around.

The staircase is the first question’s answer. It is exactly as bad as it looks and exactly as good: fifteen or sixteen decimal digits, relative, per operation, everywhere in the normal range. That is an extraordinary amount of accuracy, and the reason so much of what follows is about losing it is that losing it is easy and the ways of doing so are not obvious.

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. 6 What the second question can do to the first. A thirteen-by-thirteen system whose exact answer is the integers one to thirteen, solved in double precision — sixteen digits of relative accuracy per operation, and a result with 0.05 where 8 belongs and −2.6 where 10 belongs. The backward error of that solve is 2.2·10⁻¹⁷. Nothing in the arithmetic went wrong.

That figure is the reason this site exists, and it is the reason an essay about the spacing of floating-point numbers has to be followed immediately by one about what an error is. Sixteen digits went in. None came out. Both of those statements are true at once, and the vocabulary for saying so — forward error, backward error, condition number — is the subject of the next field.

Before that, two more essays finish the arithmetic. One is about the subtraction that reveals inherited error rather than creating it, which is where a working numerical analyst spends more time than anywhere else. The other is about the fact that the sum of a list depends on the order of the list, which sounds like a curiosity and is the difference between a simulation that conserves energy and one that does not.

Neither is about matrices. Both are load-bearing for every matrix argument that follows, because a matrix computation is nothing but a great many sums of products, performed in a particular order, and the two things this field is about are what a sum of products costs and what the order does to it.

Where this goes next

10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²110⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹xrelative error of the computed value(1 − cos x)/x², as written2 sin²(x/2)/x²no digits left at alldouble precision throughoutone function, two spellings Cancellation takes the answer, not a digit Subtracting two nearly equal numbers is exact. That is what makes it dangerous — the subtraction introduces no error at all, it exposes error the operands were already carrying, and the exposure can consume every significant figure at once. 10¹10²10³10⁴10⁵10⁶10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹number of terms addedrelative error against the exact sumin orderin a treecompensatedbinary32 · terms are 1/icompensated: 3·10⁻⁸ The order they are added in Addition is associative in the algebra and is not associative in the arithmetic. The same million numbers, added in a different order, give answers that differ in the third significant figure — and the fix is not a wider float, it is a different order. 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 The exact answer to a nearby problem A good algorithm does not give an approximate answer to your problem. It gives the exact answer to a problem very close to yours — and once that is the definition, a wrong result has two possible authors and they can be measured apart.