Numerical Linear Algebra
where the algebra stops being the arithmetic
The answer that depends on the machine

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.
33 min read 27 figures Where the algebra stops being the arithmeticMeasured, not assertedWhose fault is it

Worth reading first: Rank is a decision · The same program, twice · A block nobody can call sparse.

Rank is a decision is one of this site’s earliest essays and its claim is that rank is not a property a floating-point matrix has: the singular values decay, there is no gap, and where they are cut is a choice that should be reported with the number.

This essay is that claim with one word added. The cut is a choice, and which side of it a singular value falls on is not entirely the caller’s.

The numerical rank of one matrix is 10, 12, or somewhere betweenSingular values of one 60 × 14 matrix whose spectrum falls by a fixed factor per index, so there is no gap anywhere and the threshold decides the rank. Each curve is the same matrix, its Gram matrix formed with the inner products summed in a different number of pieces. The curves lie on top of each other for the first several values and separate below about 7.3·10⁻¹³, which is where forming AᵀA has put the rounding. The dashed line is the threshold, σ₁ · 10⁻¹⁴. Counting the values above it gives 12 at p = 1, 12 at p = 2, 12 at p = 4, 10 at p = 8, 10 at p = 16, 11 at p = 32, 11 at p = 60 — the rank of one matrix, as a function of how many workers were available. Not a digit of an answer: the number of columns a model built from this matrix would have.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
Fig. 1 The singular values of one 60 × 14 matrix, computed seven times. Each curve is the same matrix, its Gram matrix formed with the inner products summed in a different number of pieces; the dashed line is the threshold, σ₁ · 10⁻¹⁴.

Seven curves that lie on top of each other for the first several values and separate below about 10⁻¹³ of σ₁. Counting the values above the threshold gives 12, 12, 12, 10, 10, 11, 11.

One matrix. Three different answers to how many columns does this have.

What is fixed and what is not

The matrix is fixed, exactly: 60 × 14, with singular values falling by a factor of 10^1.1 per index, generated once from a seeded stream. The threshold is fixed: σ₁ · 10⁻¹⁴, which is a conventional choice and the one this site uses elsewhere. The precision is fixed at binary64.

What varies is how the inner products building the Gram matrix AᵀA were summed — one accumulator, two, four, and so on to sixty. That is the only difference between the seven curves, and it is a difference no part of the program specifies.

Why the Gram matrix and not the SVD directly

A fair objection: a rank computation would take the singular value decomposition of A rather than the eigenvalues of AᵀA, and forming the Gram matrix squares the condition number, which the road that squares the problem spends a whole essay warning against.

Both halves of the answer matter.

It is what a great deal of software does. A rank test inside a covariance routine, a least-squares solver written on the normal equations, a kernel method, a principal-component computation on a correlation matrix: all of them start from a Gram matrix, because the Gram matrix is what the problem hands them.

And the squaring is what makes the effect visible at a threshold rather than at the last bit. Squaring the conditioning halves the number of digits, so a disagreement that would sit at 10⁻¹⁶ in the singular values of A sits at 10⁻⁸ in the eigenvalues of AᵀA and therefore at 10⁻⁸ in their square roots. The effect is real in either route; the Gram route brings it up to where a conventional threshold can see it.

The honest statement is therefore: this figure measures a rank test on a Gram matrix, that is a common and slightly disreputable thing to do, and the disreputable part is what makes the machine dependence visible rather than what creates it.

The condition, stated exactly

The disagreement is not universal and the boundary is sharp enough to state.

The curves separate below the level at which the computed singular values are made of rounding. Above it every partitioning agrees to many digits. So the rank is machine-dependent exactly when the threshold falls in the part of the spectrum that is rounding — and that is the same condition under which the rank was never well defined in the first place.

The drag makes it concrete. At 0.6 decades an index the spectrum falls slowly, every singular value clears the threshold, and all seven partitionings return 14 — the full rank, with nothing to disagree about. At 1.1 the spectrum reaches the rounding level around the tenth value and the verdicts split three ways. At 1.4 and 1.6 they split further.

Which yields a rule that is more useful than the alarm: if two partitionings disagree about the rank, the threshold is in the noise, and the rank was a coin flip before the machine was involved. The machine dependence is a symptom of an ill-posed cut rather than a cause of a wrong answer.

Three thresholds, and which of them survive

The threshold in the figure is σ₁·10⁻¹⁴, which is one of three conventions in common use. It is worth asking which of them this affects, because the answer is all of them, differently.

A relative threshold, σ₁·ε. What the figure uses. Its position depends on σ₁, which every partitioning agrees about to fourteen digits, so the threshold itself is stable and only the values being compared with it move. This is the best-behaved of the three.

An absolute threshold. Used where the data has a known noise level — the discrepancy principle’s shape, which when the answer is a choice argues is the only principled way to pick one. Equally affected, and with the same structure: the threshold is fixed and the values wander across it.

A gap-based rule — cut where σᵣ/σᵣ₊₁ is largest. This one is affected worse, and it is the rule that sounds most principled. Below the noise floor the computed singular values are rounding, so the ratios between them are ratios of rounding, and the largest of them lands wherever it lands. A gap rule applied to a spectrum with no gap does not return an arbitrary-but-stable answer; it returns the position of the largest accidental ratio in the noise, which is as machine-dependent as anything in this field.

The ordering is worth carrying: a fixed threshold on a moving value is a coin flip near the line, and a rule that searches for a feature in the noise is a coin flip everywhere in it.

Why that is still worth knowing

A reader might take the previous paragraph as a reason to stop worrying. Three reasons not to.

The verdict is consumed as an integer. A rank of 10 and a rank of 12 do not produce answers that differ slightly; they produce models with a different number of columns, projections onto subspaces of different dimension, and — in a regularised solve — a truncation at a different index. The downstream difference is not proportional to the ambiguity upstream.

Nothing reports the ambiguity. A routine that returns 11 returns 11. It does not return 11, and the tenth and twelfth singular values are within the noise of the threshold, which is the information the caller needs, and which the cheap rank and what it cannot see shows is already missing from the cheap rank-revealing routines for a different reason.

And the ill-posedness is the normal case. A matrix with a clean gap in its spectrum is a matrix somebody constructed. Real spectra decay smoothly — spectral decay measures it on kernel matrices, and the reduction field’s Hankel singular values do the same thing — so the threshold usually does land somewhere in a continuum, and the question of which side a value falls on is usually being decided by something.

Where the curves separate, and what that level is

The seven curves are one curve for the first several values and seven for the rest, and the level at which they part is a computable quantity rather than an accident.

An entry of the Gram matrix is Σₖ AₖᵢAₖⱼ, a sum of sixty products. For the leading columns those products are all of similar size and the same sign, so the summation condition number is near one and every partitioning agrees to the last bit. For the trailing columns — the ones carrying the small singular values — the products cancel almost completely, κ climbs, and the partitionings disagree by κ·u.

The eigenvalues of the Gram matrix inherit that: a perturbation of size ε in the entries moves an eigenvalue by about ε‖A‖², so the computed σ² are uncertain by about u‖A‖² and the computed σ are uncertain below about σ₁√u ≈ σ₁ · 1.5·10⁻⁸.

That is the floor the square-root route imposes and this site has already measured it in another field: the product nobody had to form finds Hankel singular values flooring at σ₁√u — predicted 2.3·10⁻⁹, measured 9.7·10⁻¹⁰ — for exactly this reason. The two essays are the same arithmetic from two directions: there the loss is an accuracy problem, here it is a reproducibility one, and both come from forming a product whose conditioning is the square of the original’s.

Which locates the disagreement precisely. Below σ₁√u the computed singular values are noise, and noise is what differs between partitionings. A threshold below σ₁√u on a Gram-based rank test is reading the machine. The conventional threshold σ₁·10⁻¹⁴ is six orders below that line, which is why this figure has anything in it at all.

What to report instead

The same recommendation the field keeps producing, in its rank-shaped form.

Return the spectrum, not the count. A caller who has the singular values can apply any threshold and can see how many values are near it. This is what rank is a decision argues for on grounds that have nothing to do with machines, and the machine dependence is a second argument for the same thing.

Return the gap ratio at the cut. σᵣ/σᵣ₊₁ is one division and it says whether the cut is a cut or a guess. A ratio of 10³ is a decision; a ratio of 1.2 is a coin flip, and a coin flip is exactly where two partitionings will disagree.

And if a single integer must be returned, make the reduction order-independent. Pre-rounded accumulation on the Gram entries costs three operations an element — for a 60 × 14 matrix, nothing at all — and makes the returned integer a function of the matrix and the threshold, which is what a caller assumed it was. It does not make the integer right, because there is no right one; it makes it stable, which is what a caller who is going to store it in a model actually needs.

What a downstream model does with it

The reason a rank is worse to have unstable than a residual is that a rank is not consumed as a number. It is consumed as a shape, and shapes propagate.

Take the ordinary chain: compute a rank r, keep the leading r singular vectors, project the problem onto them, solve the reduced problem, project back. Two machines with r = 10 and r = 12 do not produce answers 20% apart — they produce answers of different quality on different subspaces, and the difference depends entirely on how much the eleventh and twelfth directions matter for the right-hand side in question.

If those directions are noise, both models are equally good, and the extra two columns in one of them are fitting rounding — which is the whole subject of the regularisation field, and when the answer is a choice is the essay about what happens when a solution is determined by a decision rather than by data.

If those directions carry signal, the smaller model is worse, and by a margin that is not small.

The uncomfortable part is that the machine has no opinion about which. It cuts where the arithmetic happened to land, and whether that was the right place is a question about the problem that the rank test was supposed to be answering.

So the composition is: an ill-posed decision, made by rounding, consumed as a shape, with downstream consequences of unpredictable size and no record anywhere that a decision was made. That is the worst shape in this field, and it is the reason the rank essay sits at the end of the verdict sequence rather than at its start.

The other quantity with the same shape

Worth naming, because it is the same measurement in the field next door.

A rank that is a number of digits makes the point that a numerical rank at a threshold ε is a statement about how many digits the data has. That reading survives here intact, and this essay adds a term to it: the rank is a statement about how many digits the data has and about how many digits the arithmetic delivered, and the second is a property of the machine.

The two are not independent. A matrix whose data carries eight digits and whose Gram matrix is computed to eight digits has a rank that is a coin flip at the eighth; the same matrix with an order-independent Gram computation has a rank that is a stable but arbitrary integer at the eighth. Neither is right. The second is at least the same everywhere, which is the property a pipeline needs even when the number in it is a convention.

One line

A rank is an integer produced by comparing a computed number with a threshold, and when the threshold sits where the arithmetic ran out, the integer is a property of how the work was divided — which is a fact about the machine appearing in the shape of a model.

What the figure would look like on the honest route

A last note to keep the essay from over-claiming, and it is a prediction rather than a measurement, which is why it is labelled as one.

Take the singular values of A directly, by a one-sided Jacobi sweep rather than through AᵀA. Nothing has been squared, so the level at which the computed values are made of rounding is σ₁·u rather than σ₁√u — eight orders lower. A threshold at σ₁·10⁻¹⁴ then sits above that level rather than below it, and the argument this essay makes has nothing to bite on: the values being compared with the threshold are determined to many digits, and a disagreement in their last bits does not move an integer.

That is a prediction from the same mechanism rather than a figure, because the site’s Jacobi sweep does not take a reduction as an argument and instrumenting it would be a second experiment. What is measured, elsewhere, is the floor itself: the product nobody had to form finds the square-root route reaching 10⁻¹⁶ where the product route stops at σ₁√u.

So the honest scope of this essay’s alarm: it is about rank tests that start from a product, which are common and which the least-squares field already warns about for a different reason. A rank test on the matrix itself is a different computation, and the argument here does not reach it.

At other settings

The numerical rank of one matrix is 12, 13, or somewhere betweenSingular values of one 60 × 14 matrix whose spectrum falls by a fixed factor per index, so there is no gap anywhere and the threshold decides the rank. Each curve is the same matrix, its Gram matrix formed with the inner products summed in a different number of pieces. The curves lie on top of each other for the first several values and separate below about 7.3·10⁻¹³, which is where forming AᵀA has put the rounding. The dashed line is the threshold, σ₁ · 10⁻¹⁴. Counting the values above it gives 12 at p = 1, 12 at p = 2, 12 at p = 4, 13 at p = 8, 13 at p = 16, 12 at p = 32, 12 at p = 60 — the rank of one matrix, as a function of how many workers were available. Not a digit of an answer: the number of columns a model built from this matrix would have.13579111310⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹index ksingular valueσ₁ · 10⁻¹⁴, the thresholdone matrix, three rankspartitionings7lowest rank12highest rank13threshold7.3·10⁻¹⁴σ₁7.3the curves separate in the noiseand the threshold is drawn through it
Fig. 2 A slower decay, where the threshold sits higher in the spectrum and fewer values are in the noise.
The numerical rank of one matrix is 9, 11, or somewhere betweenSingular values of one 60 × 14 matrix whose spectrum falls by a fixed factor per index, so there is no gap anywhere and the threshold decides the rank. Each curve is the same matrix, its Gram matrix formed with the inner products summed in a different number of pieces. The curves lie on top of each other for the first several values and separate below about 7.3·10⁻¹³, which is where forming AᵀA has put the rounding. The dashed line is the threshold, σ₁ · 10⁻¹⁴. Counting the values above it gives 10 at p = 1, 9 at p = 2, 9 at p = 4, 11 at p = 8, 10 at p = 16, 10 at p = 32, 10 at p = 60 — the rank of one matrix, as a function of how many workers were available. Not a digit of an answer: the number of columns a model built from this matrix would have.13579111310⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹index ksingular valueσ₁ · 10⁻¹⁴, the thresholdone matrix, three rankspartitionings7lowest rank9highest rank11threshold7.3·10⁻¹⁴σ₁7.3the curves separate in the noiseand the threshold is drawn through it
Fig. 3 And a faster one, where more of the spectrum is below the level the arithmetic delivered.
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. 4 The decision this essay is about, in the essay that first argued it was one.
A randomised rank-k solve, 4 seeds a rankRelative error against the rank kept, on a logarithmic vertical axis, with a vertical bar at each rank spanning the seeds. The band is 1.68 wide at rank 8, where the method is at its worst, and 1.004 wide at rank 20, where it is at its best. The same computation on the same data returns a different answer each time.0481216202428321rank keptrelative errormedianthe answer movesspread at rank 81.7spread at rank 201best median error0.15widest where the method is worstand the bound does not say so
Fig. 5 The band a rank sits in when the data has noise in it, which is the other reason the cut is a choice.
Least squares by QR and by the normal equations in binary32Relative error of the computed coefficients against epsilon, on log axes. The QR route is a flat line near the bottom; the normal-equations route climbs and then stops, at the epsilon where the cross-product matrix becomes exactly singular.10⁻⁸10⁻⁷10⁻⁶10⁻⁵10⁻⁴10⁻³10⁻²10⁻¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹10¹ε in Läuchli's matrix (smaller ε, larger κ)relative error in the coefficientsAᵀA exactly singularnormal equationsQRκ from 1.7·10⁸ to 17the cliff is at √u = 2.4·10⁻⁴
Fig. 6 The squaring that brings the disagreement up to where a threshold can see it.
What a rank-10 approximation can achieve, by spectrumA semi-logarithmic plot of singular value against index for three spectra — geometric decay, algebraic decay, and flat — with the rank-ten approximation error marked on each.1112131415110⁻³10⁻²10⁻¹1index jσⱼσ11geometric, 0.85ʲalgebraic, j^−1.00flatmeasured, and equal to σ₁₁rank-10 error, geometric0.2rank-10 error, algebraic0.091rank-10 error, flat160×60, spectrum chosen rather than the entriesthe matrix decides, not the method
Fig. 7 A spectrum with no gap in it, which is the normal case.
What a perturbation of 10⁻⁶ moves, against the gap it is applied acrossThree curves against the gap between two eigenvalues, both axes logarithmic. Two are flat; the third rises steeply as the gap closes and then levels off.10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻⁸10⁻⁵10⁻²gap between the two eigenvalueshow far it movedthe eigenvectorsthe eigenvaluestheir plane‖E‖ / gapone perturbation, three answerseigenvalue shift, spread over the sweep1plane angle, spread over the sweep1eigenvector angle, spread1.6·10⁵the dashed line is Davis–Kahan's ‖E‖/gaptwo of the three never noticed
Fig. 8 And what happens to a decision as the gap it is reading closes.
The same 4096 numbers, added up 26 ways, 21 different answersEach dot is one reduction of one vector of 4096 binary64 numbers whose summation condition number is 1.01·10⁸: the vector is cut into p contiguous pieces, each piece summed left to right, and the pieces combined in a tree. The vertical position is the distance from the exactly rounded sum, relative. Filled dots are the static split a fixed thread count gives; open dots are unequal splits, which is what a work-stealing scheduler produces at the same p. 21 of the 26 runs returned distinct values, spanning 1.55·10⁻¹³ — 2.3·10⁷ ulps of the answer. The dashed line is the classical bound γ₍ₙ₋₁₎Σ|xᵢ|, which every one of them satisfies with 2.58·10⁴ to spare, and which is the same number for all of them: it contains n and Σ|xᵢ| and nothing about the order.110¹10⁻¹¹10⁻⁸10⁻⁵pieces the vector was divided intodistance from the exact sum, relativethe published boundκ · uone vector, one algorithmdistinct answers21runs26spread, in ulps2.3·10⁷κ of the sum10⁸bound ÷ worst error2.6·10⁴nobody chose pand no answer is the answer
Fig. 9 The disagreement in the inner products underneath all of it.
The summation condition number of six inner products this site already computesκ = Σ|xᵢ| ÷ |Σxᵢ| for six reductions taken from computations elsewhere on this site, at 64 terms: an orthogonality check between two columns of a Q, one component of a residual b − Ax on an ill-conditioned system, a nearly degenerate 2 × 2 determinant, a Hutchinson trace probe zᵀAz, a residual norm rᵀr, and a curvature pᵀAp. The range is 1 to 1.01·10¹⁷. The two at the safe end are sums of squares, where no term can cancel another; the three above 10¹⁰ are the orthogonality check, the residual and the determinant — which is to say every quantity a stopping test, an orthogonality test or an orientation test is written in. The reductions that are reproducible are the ones nobody makes a decision from.κ of the suma component of b − Ax1.01·10¹⁷ad − bc, near-degenerate3.6·10¹⁶qᵢᵀqⱼ, an orthogonality check7.39·10¹⁵pᵀAp, a curvature409zᵀAz, a trace probe41.8rᵀr, a residual norm1measured, not assumedhighest10¹⁷lowest1above 10¹⁰3terms64sums of squares are safeand nobody decides anything from one
Fig. 10 Where a Gram entry sits among the site’s reductions.
Five summation policies, 400 permutations each, and how many answers each returnedThe same 1,024 numbers, summation condition number 10⁸, presented in 400 different orders. The bar is the number of distinct binary64 values the policy returned; the figure beside it is the worst relative error it made. A single accumulator returns 303 values and eight pieces return 72; Kahan's compensated loop is the most accurate of the three at 4.96·10⁻¹⁰ and still returns 119. Pre-rounded summation returns one value, at an error of 4.91·10⁻⁶ — four orders worse than the compensated loop and the same on every machine. Exact accumulation returns one value and the right one. The comparison is bitwise because the claim is bitwise: a tolerance here would pass everything.distinct answersone accumulator303eight pieces72compensated119pre-rounded1exact1worst error 1.17·10⁻⁸worst error 4.61·10⁻⁹worst error 4.96·10⁻¹⁰worst error 4.91·10⁻⁶worst error 0bitwise, or not at allpermutations400one accumulator303pre-rounded1its error4.9·10⁻⁶compensated error5·10⁻¹⁰accuracy and agreement are different propertiesand the accurate one is not the agreed one
Fig. 11 The policy that would make the count stable.
Six ways to add up a vector, priced against what they returnBinary64 operations per element on the horizontal axis and relative error on the vertical, for one vector of 2,048 numbers with a summation condition number of 1.05·10⁶. Filled markers are the policies whose answer does not depend on the order the terms arrive in; open markers are the ones whose answer does. The cheapest order-independent policy costs 3 operations an element — a divide, a round and an add — and is 2.29·10⁻⁷ wrong; the most accurate policy here is pre-rounded, two passes at 7. Nothing sits in the bottom left. The cost of an answer that is the same on every machine is between three and twelve times the cost of an answer that is not, and the two-pass pre-rounded policy is the one worth knowing about: seven operations, order-independent, and accurate.02468101210⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸binary64 operations per elementrelative errorruntime orderfixed partitionpre-rounded, one passpre-rounded, two passescompensated (Kahan)exact accumulationwhat determinism costsreproducible policies3cheapest, ops/element3its error2.3·10⁻⁷two passes10⁻¹⁷ordinary, 8 pieces2.7·10⁻¹¹nothing is cheap and exactthe middle of the figure is the answer
Fig. 12 And what it costs on a matrix this size, which is nothing.
One system, one tolerance, and 11 different amounts of workConjugate gradients on a 200 × 200 symmetric positive definite matrix with κ = 10⁴, stopped when the relative residual falls below 10⁻¹⁰, run at each of 13 partition counts. The only thing that differs between the runs is how many pieces the two global inner products were summed in — which on a real machine is the number of workers. Iteration counts run from 674 to 690, 11 of them distinct, a spread of 2.4% of the work. Every run converged and every answer is right to the accuracy asked for; the forward errors span 6.31·10⁻¹¹ to 1.08·10⁻¹⁰, a factor of 1.72. What differs is not a digit at the end of an answer — it is the bill, and it is an integer.110¹670673676679682685688691694pieces the inner products were summed initerations to the tolerance674, the cheapest run690, the dearestthe same solve, pricedpartitionings run13distinct counts11spread, per cent2.4best forward error6.3·10⁻¹¹worst1.1·10⁻¹⁰one matrix, one toleranceand the cost is the machine's
Fig. 13 The same mechanism deciding a step count.
Two conforming builds, one matrix, and 26 different verdicts in 200 triesForty 6 × 6 Gram matrices at each of five conditionings, each formed exactly from a tall matrix whose last column is nearly the sum of its first two. Every one of them is asked *are you positive definite* by a Cholesky, once with the inner products accumulated with a single rounding per term and once with each product rounded before it is subtracted. The bar is how many of the forty came back with different verdicts: 7 at 10⁻⁸, 2 at 10⁻⁹, 4 at 10⁻¹⁰, 10 at 10⁻¹¹, 3 at 10⁻¹². The two shaded parts split those by which build was right, settled by the signs of the leading principal minors computed in BigInt rationals — 9 to the fused form and 17 to the rounded one. So neither is the correct implementation. The verdict is a property of the build, and the answer to which build is right is a property of the matrix.-12-11-10-9-8024681012log₁₀ of how nearly dependent the columns areverdicts that disagreed, of 40724103which one is correctmatrices tested200verdicts disagreed26fused was right9unfused was right17lower part: thefused build was righta sign has no last digitso a verdict has nowhere to hide
Fig. 14 And a sign.
How many columns a decade of accuracy costs, measured and predicted, at q = 0.500The number of singular values above ε, against the number of digits ε asks for, on a 128 × 128 block between two intervals separated by a gap of 1. The measured curve is a straight line at 0.55 columns a decade: a digit costs the same handful of columns wherever you buy it, which is why the accuracy is a knob and not a cliff. The upper line is what the geometry alone promises — q^(p+1)/(1 − q) below ε, solved for p, using four numbers and no entry of the matrix — at 3.32 columns a decade. Both are straight and they are not the same straight: the bound is right about the shape and loose by 5.3× about the constant, which is the safe direction for a quantity you have to allocate storage from.0369121501122334455digits asked for, −log₁₀ εcolumns keptwhat the geometry promiseswhat the matrix costsa rank is a number of digitscolumns a decade0.55bound, a decade3.3rank at 10⁻⁸5bound at 10⁻⁸28q0.5the shape is rightand the constant is not
Fig. 15 What a rank verdict is worth when the gap is wide, which is when nothing here applies.
Error of the best rank-k approximation to a 12×12 matrixApproximation error against k on a logarithmic axis, with the measured error and the next singular value drawn as separate curves lying exactly on top of one another.123456789101110⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²1rank k of the approximation‖A − Aₖ‖measured, 2-normσₖ₊₁, from theorymeasured, Frobeniusthe first two agreeto 2·10⁻⁸worst |‖A−Aₖ‖₂ − σₖ₊₁| / σₖ₊₁2·10⁻⁸worst Frobenius discrepancy2·10⁻⁸κ = 10⁹; 30 random rank-3 matrices, none closerthe error is σₖ₊₁
Fig. 16 The approximation a rank decision is usually made in service of.
The rank of an admissible block and of a touching one, against how finely they are sampledBoth blocks are of the kernel 1/r; both are 32, 64, 128, 256 points a side; both are truncated at 10⁻⁸. The admissible pair — [0, 1] against [2, 3] — needs 5, 5, 5, 5 columns, which is one number. The touching pair — [0, 1] against [1, 2] — needs 9, 11, 12, 13, climbing by about one per doubling, which is a logarithm. Neither of them grows like the block, and only one of them stops. That difference is what the admissibility test in a partition is buying, and it is why the touching pair is kept dense rather than compressed at all.45678903691215log₂ of the points a sidecolumns above 10⁻⁸two intervals that touchtwo that do notthe measurement that is a null resultadmissible, n = 325admissible, n = 2565touching, n = 329touching, n = 25613stored ⁄ dense at largest0.039the rank belongs to the geometryand not to the sampling
Fig. 17 The storage that follows from the integer this essay says is unstable.
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. 18 The distinction that survives all of this untouched.
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. 19 And the standing that would settle it, at a cost nobody pays for a rank.
The smallest defect a regression test can catch, on a machine that will not repeat itselfA 150 × 150 system with κ = 10⁴, solved by conjugate gradients to 10⁻¹⁰ on six partitionings. The shaded band is how far the answer moves across those six when nothing at all is wrong: 3.24·10⁻¹², relative. Each pair of dots is one defect — a single entry of the matrix changed by a relative δ, an assembly slip or a stale coefficient — solved on the same six machines, showing its nearest and furthest run. A defect is catchable only when its nearest run clears the band, because a test whose threshold sits inside the band fails a correct build. Everything at or below 10⁻¹³ overlaps and is invisible; 10⁻¹² is the smallest that separates, at 4.61·10⁻¹². So the tolerance exists, it is bracketed between 3.24·10⁻¹² and 4.61·10⁻¹² — a factor of 1.42 — and it is neither zero nor the 10⁻⁸ that usually gets typed.-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
Fig. 20 What a test suite can ask of a verdict that moves.
What the runs disagree by, and what the bound says they mightOver ten decades of summation condition number, on vectors of 2048 numbers: the measured spread between seven partitioned reductions of the same vector, the quantity κu, and the classical error bound γ₍ₙ₋₁₎Σ|xᵢ| relative to the sum. The measured spread and κu are one curve times a constant — the ratio runs from 0.253 to 0.258 across the whole sweep, a spread of 1.021 while both quantities move by ten orders. The bound sits 7932 to 8100 above the disagreement and is identical for every ordering, which is what makes it useless for the question a reader actually has: not how wrong is this, but will the next run say the same.10²10⁴10⁶10⁸10¹⁰10¹²10⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹κ of the sumrelative sizethe boundκ · umeasured spreadtwo curves and one constantspread ÷ κu, low0.25spread ÷ κu, high0.26bound ÷ spread, low7932decades swept10the spread is computablethe bound cannot see the order
Fig. 21 The size of the disagreement being amplified into an integer.
Loss of orthogonality against condition number, in binary64A log–log plot of the norm of Q-transpose-Q minus the identity against condition number. Classical Gram–Schmidt rises steeply, modified Gram–Schmidt rises gently, and Householder is flat.110²10⁴10⁶10⁸10¹⁰10¹²10⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹condition number κ(A)‖QᵀQ − I‖classicalmodifiedHouseholderκ²uκu8×8, eight seeds per κ, binary64all three reconstruct A
Fig. 22 Orthogonality as a number, computed from an inner product at the wrong end of the census.
QᵀQ from classical Gram–Schmidt and from Householder on the 8×8 Hilbert matrixTwo eight-by-eight tables of QᵀQ. The upper one has ones on the diagonal and entries as large as one off it; the lower one is the identity to three decimal places everywhere.A = H8 · κ = 1.5·10¹⁰ · both factorisations reconstruct A to 6·10⁻¹⁷the diagonal is 1 in both — every column is a unit vector either way1.000000000001.000000000001.000000000001.000000000001.00000.002-0.002000001.0000.125-0.13300000.0020.1251.000-1.0000000-0.002-0.133-1.0001.000classical Gram–Schmidt1.000000000001.000000000001.000000000001.000000000001.000000000001.000000000001.000000000001.000Householderclassical ‖QᵀQ − I‖1.4Householder ‖QᵀQ − I‖1.4·10⁻¹⁵largest off-diagonal 1 against 3.1·10⁻¹⁶length is not angle
Fig. 23 The same quantity entry by entry, where the cancellation that makes it fragile is visible.
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. 24 A bound whose worst case is never attained, in the field that measured it first.
The backward error of a hierarchical solve, against the error of the representation it usedSolve A_H x = b exactly and the residual against the matrix that was wanted is b − Ax = (A_H − A)x, so ‖b − Ax‖ ⁄ ‖A‖‖x‖ cannot be anything but the representation's own error. Measured across six accuracies spanning ten decades, the two track at a slope of 0.994 and sit a constant 6.7× apart, which is the difference between a norm of a matrix and a norm of that matrix applied to one vector. The reading is the one this field is built on: the accuracy knob is not an accuracy, it is a backward error chosen in advance. Everywhere else on this site a backward error is something an algorithm produced and somebody then measured; here it is a line in the program, and its size is known before the solve starts.10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³‖A − A_H‖ ⁄ ‖A‖, the representation‖b − Ax‖ ⁄ ‖A‖‖x‖, the solveequala backward error, chosenslope0.99representation at 10⁻⁸1.2·10⁻⁹backward error there1.3·10⁻¹⁰representation at 10⁻¹²1.2·10⁻¹³backward error there2.5·10⁻¹⁴the compression is not an approximationit is a perturbation of the problem
Fig. 25 A backward error chosen rather than incurred, which is the reading this field gives the machine’s own.
Two perturbation bounds and the error that was measured, on a 8×8 matrix spread over 8 decades of unitsThree curves against the size of an entrywise relative perturbation. The normwise bound κ_∞·ε is a valid bound and sits 4·10⁷ times above the componentwise one cond(A, x)·ε, which is also a bound and is nearly attained by the worst of forty random perturbations at each size.10⁻¹⁴10⁻¹³10⁻¹²10⁻¹¹10⁻¹⁰10⁻⁹10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1relative size of the entrywise perturbationrelative forward errorκ_∞ · εcond(A,x) · εmeasuredboth bounds holdκ_∞(A)1.9·10⁸cond(A, x)4.8ratio of the bounds4·10⁷both curves above the data are boundsand only one of them is a measurement
Fig. 26 An error measured entry by entry rather than in a norm, which moves the quantity to the dangerous end of the census.
Backward error with and without the structure, for Levinson and for elimination, n = 12Both solvers are backward stable on this family by the usual measure: the smallest perturbation of any kind that explains either answer stays below 3.1·10⁻¹⁷ at every ρ, and the two curves are indistinguishable. Insisting the perturbation be a symmetric Toeplitz matrix — the kind of object the problem was posed with — lifts both by orders, and at ρ = 0.999 Levinson's rises by a factor of 2.48·10⁶. Neither number is alarming in absolute terms; what is worth carrying is that the reassuring one is the one that is reported, and the two are not measuring the same thing.10⁻³10⁻²10⁻¹110⁻¹⁹10⁻¹⁷10⁻¹⁵10⁻¹³10⁻¹¹10⁻⁹1 − ρbackward errorLevinson, Toeplitzelimination, ToeplitzLevinson, any kindelimination, any kindat ρ = 0.999Levinson, any perturbation3.1·10⁻¹⁷Levinson, Toeplitz only7.6·10⁻¹¹the ratio between them2.5·10⁶diagonal defect of the first0.97the number that is reportedand the number that is asked about
Fig. 27 A backward error measured against a structure rather than against any one matrix.

Where this goes next

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 Rank is a decision A floating-point matrix does not have a rank. It has a spectrum of singular values, and somewhere in that spectrum is a place where the values stop being signal and start being noise. Deciding where is a judgement, and the evidence for it is a gap. 110¹670673676679682685688691694pieces the inner products were summed initerations to the tolerance674, the cheapest run690, the dearestthe same solve, pricedpartitionings run13distinct counts11spread, per cent2.4best forward error6.3·10⁻¹¹worst1.1·10⁻¹⁰one matrix, one toleranceand the cost is the machine's A stopping test is a race One matrix, one right-hand side, one tolerance, thirteen partition counts — and eleven different iteration counts between 674 and 690. Every run converged, every answer is right to the accuracy asked for, and what differs is the bill. 171319253137434910⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²110²indexmagnitude|r_kk|σ_kone factorisation, two verdicts‖AP − QR‖/‖A‖1.1·10⁻¹⁶|r_nn|8.7·10⁻⁴σ_min3.7·10⁻¹²column interchanges0|r_nn| is never below σ_minso the cheap verdict errs one way only 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.

What links here

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

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.

Condition squaringGram matrixNumerical rankRankReduction orderRun-to-run variationSingular valuesSpectral decay