The inverse that is never formed

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.

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¹⁴:

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

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.

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ⱼ)xⱼ = eⱼ for some small ΔAⱼ — 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.

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.

And the fingerprint sharpens with the conditioning rather than staying put.

Only one of the two sides moves. Across n = 6, 8, 10, 12 and 14 the left residual runs 8.1·10⁻¹⁸, 4.5·10⁻¹⁸, 5.1·10⁻¹⁸, 3.5·10⁻¹⁸ and 3.3·10⁻¹⁸ — flat, and if anything drifting down — while the right one runs 5.2·10⁻¹⁷, 3.7·10⁻¹⁷, 4.2·10⁻¹⁶, 6.8·10⁻¹⁶ and 1.3·10⁻¹⁵. The ratio goes 6, 8, 82, 195, 397.

So the direction the computation ran in is recoverable from the output, and increasingly so. At κ = 10⁷ a caller comparing the two residuals would see a factor of six and call it noise; at κ = 10¹⁷ they would see a factor of four hundred. The quantity that ought to be zero on both sides by definition is instead a record of the loop order, legible in the answer, and it is the cheapest available proof that X̂ is not the inverse of anything — two matrix products, no factorisation, and a number that grows precisely when it matters. A library that returned both would be telling its caller something no documented error bound does.

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 obvious repair is to check the other side instead — the one nobody computes, because it costs another n³. That repair does not work either, and the next section is why.

The expensive diagnostic does not work either #

Sweep the condition number and compute both residuals properly — each scaled by ‖A‖‖X̂‖, as a residual has to be — on 30 × 30 matrices:

at κ = 10² they are 3.3·10⁻¹⁷ and 9.1·10⁻¹⁷. At 10⁸, 1.4·10⁻¹⁷ and 1.2·10⁻¹⁶. At 10¹⁴, 1.3·10⁻¹⁷ and 7.6·10⁻¹⁶. The ratio between them runs 2.8, 9.0, 59.0 — so the asymmetry is real, it grows with the conditioning, and it is a factor of tens.

Both sides sit at the unit roundoff at every condition number. And the instability they are supposed to reveal, at κ = 10¹⁴, is 4.5·10⁻⁵.

Eleven orders of magnitude. The expensive residual — the one that costs a second n³, the one this essay has just called the side carrying information — comes back at 7.6·10⁻¹⁶ on a matrix where solving through the inverse has a backward error of 4.5·10⁻⁵. It is not a weak signal. It is no signal.

The reason is the same construction argument run once more, and it applies to both sides. ‖X̂A − I‖ measures whether X̂ is a left inverse of A, and X̂ is very nearly one — each of its columns is the exact solution of a nearby system, and nearby systems have nearly the same inverse. What is wrong with X̂ is not that its entries are far from A⁻¹’s in any norm that ‖A‖‖X̂‖ does not already divide out. What is wrong is that its columns are the inverses of different nearby matrices, and no residual of X̂ against A can see that, because a residual of X̂ against A is a statement about X̂ as a whole.

The quantity that does see it costs 2n² #

It is the backward error of the solve: ‖b − Ax̂‖ divided by ‖A‖‖x̂‖ + ‖b‖. Across the same sweep it returns 2.4·10⁻¹⁶, 2.4·10⁻¹¹ and 4.5·10⁻⁵ against the LU route’s flat 10⁻¹⁷, which is the measurement this essay opens with.

One matrix–vector product and two norms. Cheaper than the multiply it is checking, let alone than the n³ of the residual that does not work.

So the honest form of this section is not that a code checks the wrong side of an identity. It is that forming X̂ turns a solve into a multiply, and nobody computes the residual of a multiply. A programmer who has just written x = Xinv @ b is not in the frame of mind that asks what system does this x solve; the question does not arise, because there is no visible system. That is a property of the notation rather than of the arithmetic, and it is the mechanism by which a defect measured at 4.5·10⁻⁵ goes unreported in code that is otherwise careful.

It also says what a library could do about it, which is nothing clever. A routine returning an explicit inverse could return, beside it, the residual of one representative solve — or simply refuse to be the last step, and hand back an object that computes x and its residual together. The site’s standing habit is that a quantity deciding something should be printed beside the thing it decided, and here the quantity costs 2n² and the thing it decides is whether the answer means anything.

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.

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.

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.

One more reading of the sweep, because it settles which of the two routes the conditioning is responsible for. The LU route’s backward error across κ from 10² to 10¹⁴ is 2.2, 2.9, 3.3, 3.8, 2.9, 2.3 and 2.2 times 10⁻¹⁷ — flat, and if anything slightly lower at the hard end than at the easy one. The inverse route’s is 2.4·10⁻¹⁶, 2.4·10⁻¹⁴, 1.9·10⁻¹², 2.4·10⁻¹¹, 6.4·10⁻⁹, 6.8·10⁻⁷ and 4.5·10⁻⁵. One of those two sequences is a property of the algorithm and the other is a property of the matrix being multiplied through, and the flat one is the algorithm’s.

That is worth stating because the usual defence of the inverse route is that both routes lose accuracy on an ill-conditioned matrix, so the difference is academic. Both routes do lose accuracy — the forward errors are 2.8·10⁻⁴ and 1.5·10⁻², which is a factor of fifty rather than twelve orders. What separates them is not how wrong the answer is. It is whether there exists a nearby problem the answer is exactly right for, and the whole apparatus of blame rests on there being one.

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.

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.

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.

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.

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.