Numerical Linear Algebra
where the algebra stops being the arithmetic
Elimination, and the swap

The inverse that is never formed

x = A⁻¹b is how the solution of a linear system is written and it is not how it is computed. The usual reason given is cost — three times the arithmetic. The real reason is that one of the two routes is backward stable and the other is not, and at κ = 10¹⁴ they differ by twelve orders of magnitude in the number that says whose fault a wrong answer is.
21 min read 15 figures Whose fault is itIdentical algebra, different arithmeticMeasured, not assertedTwo routes to a numberAssertions that reject

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

Every numerical computing course contains the sentence do not invert the matrix, and almost every one gives the wrong reason for it.

The reason given is cost. Forming A⁻¹ takes 2n³ operations against the (2/3)n³ of an LU factorisation, so it is three times the work; and applying the result to a vector is 2n² operations, which is the same as a pair of triangular solves. Three times the arithmetic for no saving. It is true, and it is a factor of three, and a factor of three has never stopped anyone.

The real reason is that the two routes are not both backward stable, and on this site that is not a technicality — it is the sentence the whole collection is built around. A backward-stable method returns the exact answer to a nearby problem, which means that when its answer is wrong there is a way to say whose fault it is. One of these two routes gives up that right.

Backward and forward error of two routes to x, on 30×30 systems across eight decades of κBoth routes start from the same LU factorisation. The LU solve's backward error is 2.2·10⁻¹⁷ at every conditioning — a flat line at the unit roundoff. Multiplying by the explicitly formed inverse gives 4.5·10⁻⁵ at κ = 10¹⁴, a slope of 0.94 against κ. The two forward errors, drawn above them, are 2.8·10⁻⁴ and 0.015 — within a factor of 54, which is why the difference between the two methods is invisible to anyone measuring how wrong the answer is.10²10⁴10⁶10⁸10¹⁰10¹²10¹⁴10⁻¹⁸10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1κ₂(A)relative errorforward, A⁻¹bforward, A\bbackward, A⁻¹bbackward, A\bat κ = 10¹⁴η, LU solve2.2·10⁻¹⁷η, via the inverse4.5·10⁻⁵forward, LU solve2.8·10⁻⁴forward, via the inverse0.015one factorisation, two ways to use itand one of them forfeits the backward error
Fig. 1 Both routes start from the same LU factorisation. The two lower curves are the backward errors: one flat at the unit roundoff across eight decades of conditioning, one on a slope of one. The two upper curves are the forward errors, which sit almost on top of each other. Drag the size of the system.

The experiment, which controls for everything else

The comparison is arranged so that nothing but the last step differs.

A 30×30 matrix with a prescribed condition number is factorised once, by the same luFactor in both runs. Route one solves against b. Route two solves against each column of the identity, assembles the result into X̂, and multiplies X̂b. Same factorisation, same pivots, same arithmetic; the only difference is what the triangular solves are aimed at.

The right-hand side is built from a chosen x, so the forward error is known rather than estimated — the site’s standing habit, and the only reason the last two columns of the table below mean anything.

At κ = 10¹⁴:

LU solve inverse, then multiply
backward error η 2.2·10⁻¹⁷ 4.5·10⁻⁵
forward error 2.8·10⁻⁴ 1.5·10⁻²

The forward errors differ by a factor of fifty. The backward errors differ by twelve orders of magnitude.

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: κ · u30×30, 20 seeds per κ; dashed is the worstthe problem worsens, not the method
Fig. 2 The two errors, drawn as the site first drew them. A reader who measures only the vertical axis here would call the two methods much of a muchness and blame the condition number.

The slope is the mechanism

The inverse route’s backward error is not merely larger; it grows as the first power of the condition number. Fitted across the decades where the effect is above rounding, the slope is 0.94.

That number names the mechanism. κu is the forward error of the inversion — how wrong the entries of X̂ are. Multiplying by a matrix whose entries are wrong by a relative κu produces a vector whose residual is κu‖A‖‖x‖. So the second route’s backward error is the first route’s forward error, which is a sentence worth reading twice: the inversion’s inaccuracy is inherited by the multiply as instability.

The LU route’s curve is flat at 10⁻¹⁷ across the whole axis. That flatness is what backward stability looks like as a measurement: a bound that does not depend on the matrix.

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. 3 The identity the two curves sit inside. Forward error ⪅ condition number × backward error — and when one method’s backward error is κ times another’s, their forward errors differ by a constant while their claims about whose fault it is differ completely.

The computed inverse is not the inverse of anything

There is a cleaner way to say what has gone wrong, and it has a measurement of its own.

Each column of X̂ is a backward-stable solve. So column j is the exact solution of (A + ΔA_j)x_j = e_j for some small ΔA_j — a different perturbation for each column. There is no single ΔA making X̂ = (A + ΔA)⁻¹. The assembled matrix is not the inverse of any nearby matrix, and the product X̂b is therefore not the solution of any nearby system.

That is not an abstraction; it leaves a fingerprint. The true inverse satisfies both AX = Id and XA = Id equally. The computed one does not.

‖AX̂ − I‖ and ‖X̂A − I‖ for a computed Hilbert inverse, to n = 12Two residuals of the same computed inverse, each divided by ‖A‖‖X̂‖ so both are dimensionless. The exact inverse satisfies both to zero. The computed one satisfies ‖AX̂ − I‖ at 3.5·10⁻¹⁸ — the side its columns were solved along — and ‖X̂A − I‖ at 6.8·10⁻¹⁶, 195 times larger, at n = 12. Each column of X̂ is the exact solution of a slightly perturbed system, but a different perturbation for each column, so there is no single nearby matrix whose inverse X̂ is.468101210⁻²⁰10⁻¹⁹10⁻¹⁸10⁻¹⁷10⁻¹⁶10⁻¹⁵10⁻¹⁴10⁻¹³nresidual / (‖A‖ ‖X̂‖)‖X̂A − I‖‖AX̂ − I‖the two sides, comparedn = 6, ratio6.4n = 8, ratio8.3n = 10, ratio82n = 12, ratio195the exact inverse satisfies bothand the computed one satisfies the side it was computed along
Fig. 4 ‖AX̂ − I‖ and ‖X̂A − I‖ for a computed Hilbert inverse, each divided by ‖A‖‖X̂‖ so both are dimensionless. The side the columns were solved along is at the unit roundoff; the other is up to 195 times larger.

At n = 12 the left residual is 3.5·10⁻¹⁸ and the right one is 6.8·10⁻¹⁶. There is nothing in the algebra to distinguish them — AX = Id and XA = Id are the same statement about the same object. The asymmetry exists because the computation went down the columns, and the rows were never computed as anything.

It reverses on a badly scaled matrix, which is what stops this being a rule with a direction. What it is is a rule with a shape: the residual is small on the side the arithmetic went.

The asymmetry is also the reason a common diagnostic is worthless. Code that wants to know whether its inverse is any good typically computes ‖AX̂ − I‖ and finds it at the unit roundoff, which is reassuring and means almost nothing: that residual is small by construction, since each column was the output of a backward-stable solve against exactly that equation. Checking the side the computation went down is checking that the triangular solves ran. The other side is the one carrying information, and nobody computes it because it costs another n³.

Two condition numbers of one 8×8 system, as its rows are put into different unitsFour curves against the spread of the row units, in decades. κ_∞ of the scaled matrix rises from 9.83 to 1.9·10⁸ while the componentwise condition number stays at 6.98 throughout — the same system, the same solution, and one of the two numbers is a fact about the units. Hilbert's two numbers are drawn flat beside them at 3.4·10¹⁰ and 1.2·10¹⁰: a matrix whose sensitivity no scaling repairs.02468110²10⁴10⁶10⁸10¹⁰10¹²spread of the row units (decades)condition numberκ_∞(DA)cond(DA)Hilbert κ_∞Hilbert condone system, two numbersκ_∞ at no spread9.8κ_∞ at 8 decades1.9·10⁸cond, either end7Hilbert, equilibrated1.3·10¹⁰the solution is the same at every spreadand one of these curves knows it
Fig. 5 A badly scaled matrix, where the asymmetry runs the other way. The general statement survives — the two residuals differ — and the direction does not.

And there is no number of right-hand sides that pays for it

The usual defence is that the inversion is an investment: expensive once, cheap thereafter, worth it where many systems are to be solved with the same matrix.

It is arithmetically false, and the arithmetic is one line. For m right-hand sides:

factor once, apply m times :   (2/3)n³ + 2mn²
invert once, multiply m times:      2n³ + 2mn²

The m-dependence is identical. Applying stored LU factors to a vector is two triangular solves, 2n² operations; multiplying by an explicit inverse is a matrix–vector product, 2n² operations. The same count, to the flop. So the (4/3)n³ spent forming the inverse is a constant that is never recovered — the ratio falls towards one as m grows and never reaches it, and at n = 100 the gap is 1.33 megaflops at m = 1 and 1.33 megaflops at m = 1,000.

The defence is not merely wrong about the size of the saving. There is no saving of any size, at any m, ever.

Communication for a 512×12 factorisation on 16 processorsTwo counts for the same factorisation. Rounds on the critical path: 48 for the column sweep, 4 for the reduction tree, 4 for Cholesky QR — a factor of n between the first and the other two. Words sent: 1170, 1170 and 2160 — the sweep and the tree send the same number, and the method with the fewest rounds sends the most.rounds on the critical pathHouseholder sweep48reduction tree4Cholesky QR4words sentHouseholder sweep1170reduction tree1170Cholesky QR2160two counts, two rankingsrounds, sweep ÷ tree12words, Cholesky ÷ tree1.8arithmetic, tree ÷ sweep1.5the rounds separate the threeand the words do not
Fig. 6 The one setting where the flop count stops deciding, which is the honest qualification: a matrix multiply moves data differently from a triangular solve, and the cost field is about where that matters. It does not rescue the inverse, because the (4/3)n³ is still spent.
A 6.25-bit block format against 8-bit E4M3, read two waysFour curves of relative error against the spread of the data. Two, for the block format, rise steeply when read as a median and gently when read as a norm; the other two are nearly flat.048121610⁻²10⁻¹1octaves of spread within a blockrelative errorblock, medianblock, 2-normE4M3, medianE4M3, 2-normthe same data, two readingsblock format, bits a value6.3norm across 16 octaves1.8median across 16 octaves62entries deleted at the wide end383640 values, blocks of 32a norm is dominated by what a block format keeps
Fig. 7 And the crossing that governs those arguments. Nothing in it makes an object that costs three times as much to build and the same amount to use into a bargain.

What it costs to repair, which locates the defect

If the trouble were the conditioning, correcting the answer at the same precision would not help. It does.

One step of iterative refinement — compute the residual r = b − Ax̂, solve for a correction, add it — takes the inverse route’s backward error from 4.5·10⁻⁵ to 8.1·10⁻⁹. Four steps take it to 3.1·10⁻¹⁷, which is where the LU route started. Each step costs 2n², which at n = 30 is 1,800 operations against the 54,000 the inversion cost.

Iterative refinement of the inverse-and-multiply solve at κ = 10^14The backward error starts at 4.5·10⁻⁵ and falls by about κu a step — 8.1·10⁻⁹, 5.3·10⁻¹², 3.3·10⁻¹⁵, 3.1·10⁻¹⁷, 2.8·10⁻¹⁷ — reaching the LU route's 2.2·10⁻¹⁷ after 5 corrections, each costing 2n² flops against the 2n³ the inversion cost. The forward error, drawn above, does not improve: 0.015 to 7.1·10⁻⁴, against the LU route's 2.8·10⁻⁴. What refinement at working precision buys is stability, and the accuracy floor belongs to the problem.01234510⁻¹⁸10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1correction steprelative errorLU route: η = 2.2·10⁻¹⁷LU route: forward 2.8·10⁻⁴forward errorbackward errorwhat a correction buysη before refinement4.5·10⁻⁵η after four steps2.8·10⁻¹⁷forward, unchanged7.1·10⁻⁴cost of a step, flops1800the residual is repairableand the accuracy floor is the problem's
Fig. 8 The correction steps. The backward error falls by about κu a step and reaches the LU route’s number in four; the forward error, drawn above, does not move. Drag the conditioning.

And the forward error does not move at all: 1.5·10⁻² at the start, 2.5·10⁻³ after four steps, against the LU route’s 2.8·10⁻⁴. It wanders around κu and stays there, because κu is a property of the problem and refinement performed at the working precision cannot see past it.

Both halves matter. The first says the defect is real and repairable — what was lost was the stability of one multiply, not information about the answer. The second says the repair is not an accuracy story: buying the accuracy back needs a wider precision for the residual, which is a different method.

Iterative refinement from a 24-bit factorisation, κ = 10⁴A semi-logarithmic plot of forward error against refinement step. One curve falls steeply to the level of a double-precision solve; the other is nearly flat.012345610⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹refinement step‖x − x*‖ / ‖x*‖a full double-precision solveresidual in24-bitresidual indoubleone argument apartκ·u of the factorisation6·10⁻⁴double residual, final3.2·10⁻¹³same-precision, final1.3·10⁻⁴30×30, κ = 10⁴, same factors in both runsidentical cost
Fig. 9 The mixed-precision version, where the residual is computed in a wider format and the forward error does move. That is the essay this one hands the reader to.

What the entries of an inverse are genuinely for

This is not an argument that A⁻¹ should never exist. Two places on this site need its entries, and neither is a solve.

Skeel’s condition number is ‖ |A⁻¹| |A| ‖_∞ — absolute values taken entry by entry, before any norm. That needs the entries of the inverse, not its action on a vector, and a condition number scaling cannot move forms the inverse and says so.

The variance of a least-squares estimate is σ²(AᵀA)⁻¹, and the diagonal of that matrix is what a standard error is computed from. It is the answer a user reads, not an intermediate on the way to something else.

The rule is therefore not never invert. It is: the inverse is the answer to a question about the entries of A⁻¹, and it is the wrong intermediate for a question about the action of A⁻¹ on a vector. Every case where inverting is right is a case where somebody wants to look at the numbers.

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. 10 The bound that needs those entries. Computing it requires forming |A⁻¹|, which is the legitimate case — the entries are the answer rather than a step towards one.

There is a third case, and it is the one that most often gets a working programmer into trouble: a formula written down in a paper. A Kalman gain, a projection matrix, a normal-equations estimate, a Woodbury correction — all of them are written with an inverse in the middle, because that is how linear algebra is notated, and all of them are implemented as solves. The translation is mechanical and it is not automatic: inv(R) * H * P and R \ (H * P) are the same object and different computations, and the second is the one with a bound on it.

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. 11 The best-known instance. The normal equations are written with an inverse and are usually implemented with a Cholesky factorisation, and the argument against forming AᵀA at all is a further step in the same direction — do not build the object, solve the problem.
Leverage and the deleted residual for 40 observations of a 6-column fitThe upper panel is the diagonal of the hat matrix, one bar per observation, with the average p/m = 0.150 drawn through it; the leverages sum to 6.000000000, which is exactly the number of columns. The lower panel is the leave-one-out residual, computed in closed form as eᵢ/(1 − hᵢ) and, separately, by refitting the model 40 times without each observation; the two agree to 1.3·10⁻¹². The first observation carries a leverage of 0.5000 by construction, and 1 − h is the number a hyperbolic downdate takes the square root of.the diagonal of the hat matrix, hᵢ = aᵢᵀ(AᵀA)⁻¹aᵢaverage p/m = 0.15010the leave-one-out residual: eᵢ/(1 − hᵢ), and forty refitsbars: closed form · dots: refitted without that pointone number, two fieldsΣ hᵢ, exactly p6largest leverage0.5closed form against refits1.3·10⁻¹²1 − h of the first row0.5y appears in the residualand nowhere in the leverage
Fig. 12 And the statistical one: the diagonal of a hat matrix is aᵢᵀ(AᵀA)⁻¹aᵢ, a quantity every regression package reports, and there is no way to phrase it as an action on a vector.

The shape this argument keeps taking

Three essays in this phase have now made the same move, from three unrelated directions.

A determinant is replaced by an accumulated logarithm, because the object overflows and the substitution does not. Cramer’s rule is replaced by elimination, because the formula names its products in advance and the algorithm can choose. And here the inverse is replaced by a pair of triangular solves, because the assembled object is not the inverse of anything while the solves each answer a nearby question.

In every case the object the question names is not the object worth computing, and in every case the substitute is better in a way that is measurable rather than merely cheaper. That is the sentence this phase is about, and it is worth noticing that the notation is what misleads: A⁻¹b is written as a product of two things, and it is not one.

log₁₀|det H_n| by three routes, to n = 30The exact rational determinant of the Hilbert matrix falls to 10^-518 at n = 30. The accumulated logarithm of the pivots follows it to n = 13 and then departs, reaching 10^-352 — wrong by 166 orders of magnitude and still an ordinary-looking number. The product of the pivots is a separate curve that ends at n = 29, where it underflows to exactly zero.26101418222630-550-450-350-250-150-50nlog₁₀|det H_n|exact rationalΣ log|u_ii|Π u_iithree routes, one theoremexact at the largest n-518accumulated logarithm-352decades of disagreement166smallest pivot at that n1.2·10⁻¹⁷every pivot is a normal numberat every size on this axis
Fig. 13 The same move in the determinant’s case: the object leaves the format and the substitution does not. There the substitution repaired representability and not accuracy; here it repairs stability and not accuracy, which is a different half of the same distinction.
Backward error of Cramer's rule and of elimination on 11974 2×2 systems at 24 bitsThe share of systems whose normwise backward error exceeds a given multiple of the unit roundoff. Elimination never exceeds 1.2u. Cramer's rule reaches 395u on the same systems, on a matrix whose condition number is 3.39·10⁶ and which elimination solved to 0.3u. The control curve, over uniformly random 2×2 systems, stops at 4.9u — so the failure belongs to the near-parallel family and not to a choice of scale.10⁻²10⁻¹110¹10²10³10⁴10⁻⁴10⁻³10⁻²10⁻¹1backward error, in units of ushare of systems above ituCramereliminationCramer, control24-bit arithmeticworst Cramer, in u395worst elimination, in u1.2control, worst Cramer4.9κ of the worst system3.4·10⁶one derivation, two computationsand only one of them is stable
Fig. 14 And in Cramer’s rule’s case, where the same six numbers go through two formulas and only one of them answers a nearby question.

Where the notation wins anyway

The rule do not form the inverse is one of the best-known pieces of advice in numerical computing and it is routinely ignored, which is worth accounting for rather than deploring.

The notation is a product. A⁻¹b is written as two things multiplied, and a language that lets a matrix be inverted lets that expression be typed. Every array language has an inv, and in most of them it is one character shorter than the solve. NumPy’s np.linalg.inv(A) @ b reads more naturally than np.linalg.solve(A, b) to anyone whose first language was mathematics rather than LAPACK.

A formula in a paper has an inverse in it. Kalman’s gain is PHᵀ(HPHᵀ + R)⁻¹; a projection is A(AᵀA)⁻¹Aᵀ; a Mahalanobis distance is xᵀΣ⁻¹x; a Woodbury correction has three. Every one of those is implemented as a solve, and the translation is a step somebody has to take. It is mechanical and it is not automatic, and the version with the inverse in it is what the paper says.

And the pseudoinverse is worse, because there is no notation for the solve. A⁺b is the least-squares solution, and a language that offers pinv invites forming an m×n object to apply it once. The right computation is a QR factorisation, and the site’s least-squares field is about which one; but pinv(A) @ b is what gets typed, and it costs an SVD.

The honest summary is that this is a user-interface failure rather than an ignorance one. The people forming inverses know the advice. What they have in front of them is an expression that says A⁻¹b, a library function called inv, and a translation step with no compiler to perform it. Julia’s \ operator and MATLAB’s are the design that solved it, and they solved it by making the solve the short thing to type — which is the only fix that has ever worked.

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. 15 The best-known formula with an inverse in the middle, and the two ways of computing it. The difference between them is the same shape as the difference in this essay, one level up.

It is worth noticing what the two languages that solved it did not do: neither removed inv. The function is still there, it is still correct, and it is still the right answer to a question about the entries. What changed is which expression is shorter, and that is the whole of the intervention — a default rather than a prohibition, which is the shape every successful version of this advice has taken.

What is worth carrying

The cost objection is real and it is the smaller one. Three times the arithmetic; twelve orders of magnitude in the backward error. Only the second changes what a wrong answer means.

A computed inverse is not the inverse of a nearby matrix, and the fingerprint is that it satisfies its two residuals unequally — by a factor of 195 on the Hilbert matrix at n = 12, in whichever direction the arithmetic went.

And no number of right-hand sides pays for it, because applying stored factors and multiplying by an inverse cost the same 2n² each. The investment argument is not a bad trade; it is not a trade.

The check most code performs is the one that cannot fail. ‖AX̂ − I‖ is small because that is the equation each column was solved against; the informative residual is the other one, and it costs as much as the inversion did.

The next essay keeps the same question — what does a shortcut charge? — and asks it of a formula that is genuinely worth a factor of n: a correction cheaper than the problem.

Where this goes next

10⁻²10⁻¹110¹10²10³10⁴10⁻⁴10⁻³10⁻²10⁻¹1backward error, in units of ushare of systems above ituCramereliminationCramer, control24-bit arithmeticworst Cramer, in u395worst elimination, in u1.2control, worst Cramer4.9κ of the worst system3.4·10⁶one derivation, two computationsand only one of them is stable A rule that is correct and unusable Cramer's rule gives every component of the solution in closed form, in terms of determinants, and it is a theorem. On two-by-two systems whose rows are nearly parallel it returns an answer with a backward error of 458 units of roundoff where elimination returns 1.3 — on a matrix whose condition number is 32,000 and which elimination solved perfectly. 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. 012345610⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹refinement step‖x − x*‖ / ‖x*‖a full double-precision solveresidual in24-bitresidual indoubleone argument apartκ·u of the factorisation6·10⁻⁴double residual, final3.2·10⁻¹³same-precision, final1.3·10⁻⁴30×30, κ = 10⁴, same factors in both runsidentical cost Buying the accuracy back Factorise in single precision, then correct the answer using residuals computed in double, and the result is what a full double-precision solve would have given. Compute those residuals in single instead and the identical algorithm, at identical cost, recovers nothing.

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 errorCondition numberFlop countForward errorHilbert matrixIterative refinementLU factorisationMatrix inverse