An operator with no entries

Every method in this collection so far has been handed a matrix: an array of numbers it can index, compare, pivot on, scale, count the zeros of. At the sizes where linear algebra is genuinely expensive that array does not exist.

What exists is a subroutine. Give it a vector, it returns A times the vector. A finite element code assembles the product element by element and never assembles the matrix; a spectral method transforms, multiplies by a diagonal and transforms back; a Newton method differences its own residual; a kernel method evaluates a function of pairs of points. None of them can hand over an entry, and for most of them the entries would not fit if they could.

This essay is about what survives that, what does not, and about a derivative that is six orders of magnitude worse than the exact one and costs nothing.

Those three exponents are the whole of the mechanism, and they hold at every grid the figure will draw.

The best ε has moved with the grid, which is the second column of the readout and is not something a single figure can show. It sits at 10⁻⁶ for the forward quotient at sixteen, thirty-two and sixty-four points, and at 10⁻⁵ from a hundred and twenty-eight — one decade, in the direction of a larger step.

Across n = 16, 32, 64, 128 and 200 the forward truncation slope reads 1.01, 1.01, 1.01, 1.00 and 1.00, the central one 2.00, 2.00, 2.00, 1.98 and 1.84, and the cancellation slope −0.98, −0.98, −1.00, −0.99 and −1.00. First order, second order, and one over ε: three exponents from three lines of Taylor series, none of them fitted to better than a hundredth at four of the five sizes.

The one that drifts is the central truncation slope at n = 200, and the reason is in the column beside it. The optimal ε has moved from 10⁻⁴ to 10⁻³, so the truncation branch has one fewer decade to be fitted over before the cancellation branch takes it — the exponent has not changed, the range the fit sees has.

And the floor gets better as the grid gets finer, which is the opposite of what refining a grid usually does. The forward floor reads 2.94·10⁻¹⁰, 1.2·10⁻¹⁰, 1.28·10⁻¹⁰, 3.45·10⁻¹¹ and 2.24·10⁻¹¹ across the five sizes — a factor of thirteen better at two hundred points than at sixteen — while the best ε rises by a decade. Both follow from the same fact: the floor is about √u times the ratio of ‖F‖ to ‖J v‖, and refining the grid makes the residual smaller relative to its own derivative. A matrix-free Newton method is therefore more accurate on the problems it exists for.

What survives, which is more than one expects #

Every Krylov method, unchanged. Conjugate gradients, GMRES, LSQR, MINRES, Lanczos, Arnoldi — a Krylov method’s only interaction with A is the product Av. That is not an approximation to the statement; it is the whole interface. The site’s trace estimators are matrix-free, the Krylov exponential of the previous essay is matrix-free, and the Newton–Schulz iteration of the essay before it is matrix-free, and none of them was designed to be.

That last observation is worth pausing on. Three unrelated methods in this phase turned out to require nothing but products, and in each case it was a consequence of avoiding an expensive object rather than a design goal. Avoiding the object and avoiding the entries turn out to be the same discipline.

What does not #

Everything that reads an entry, and the list is longer than it looks.

Elimination, because there is nothing to compare when choosing a pivot. Cholesky, for the same reason — including the one whose failure is the answer. Equilibration and Skeel’s condition number, which need row norms. The growth factor. The sparsity pattern, and therefore every fill-reducing ordering. Every incomplete factorisation, which is defined by a pattern. Rank, which needs an SVD.

That is the whole of the site’s scaling field, its sparse-elimination field and its rank-revealing machinery, unavailable.

What is left for preconditioning is the difficulty, because a Krylov method without a preconditioner is not a solver at scale — it is a solver at the rate the conditioning allows, which is a rate the condition number predicts and is usually too slow.

The measurement, and the awkward direction of it #

On the two-dimensional model problem at 400 unknowns:

The diagonal is available: n probes with unit vectors, or more usually the analytic diagonal, which most operators can supply. It buys nothing at all here — the diagonal of this operator is constant, so scaling by it is scaling by a scalar, and a Krylov method does not notice a scalar.

The one that helps by a factor of 2.7 is defined by the sparsity pattern.

And the ratio appears to grow with the problem: 1.93 at 64 unknowns, 2.67 at 400. Two points are not enough to say that — a ratio settling towards an asymptote passes through the same two numbers — so the sweep is below, and it does grow.

How fast the tax grows, which two points cannot say #

The theory would have allowed either answer. Unpreconditioned conjugate gradients on a two-dimensional Poisson operator takes O(√κ) iterations and κ goes as h⁻², so the count goes as h⁻¹ — that is, as the number of grid points along a side. An incomplete Cholesky improves the constant, and if it improved nothing else the ratio between the two would be flat. Whether it is flat is a question about two exponents, and two data points cannot resolve two exponents.

Eleven can:

Fitted in the log, the unpreconditioned count grows like k^0.936 — h⁻¹, as the theory says, to within the noise of eleven integers — and the preconditioned one like k^0.636. The exponents differ, so the ratio grows, and it grows like k^0.30, which is n^0.15.

That exponent is the number the section above wanted and could not have. It says the tax is real and that it is slow: a factor of sixty-four in the number of unknowns moved it from 1.9 to 3.2, and carrying the same fit out to a million unknowns gives about eight. Eight is a serious number and it is not a catastrophe, which is a more useful thing to know than increases with size — a phrase that is equally true of something heading for 3.5 and of something heading for a thousand.

It also puts the diagonal row on firmer ground. Diagonal preconditioning does not merely help little here: it returns the identical integer at all eleven sizes, 19 against 19 and 139 against 139. That is not a measurement about how good a diagonal preconditioner is. It is the statement that this operator’s diagonal is a constant, so preconditioning by it multiplies the whole matrix by a scalar, and a Krylov method’s iterates are invariant under that — the same subspace, the same minimisation, the same iterate, bit for bit. On a variable-coefficient problem the row would look completely different, which is why it appears again in the list of what a matrix-free code can still use.

One more reading is worth taking off the same table, because it changes where a matrix-free code should look for its preconditioner. The IC(0) count grows like k^0.636, which is not a constant — so the incomplete factorisation is not curing the conditioning, only improving its constant and its exponent a little. A method whose count did not grow at all would be a rate that does not notice the size, and that is multigrid rather than an incomplete factorisation. Which means the entries are buying a factor here and not an asymptote, and the thing a matrix-free code actually loses by not having them is smaller than the thing it could have had by building a geometric hierarchy instead. The list further down this essay puts multigrid second for that reason, and the exponent is why.

assertTheTaxGrowsLikeAPowerOfTheGrid holds the two exponents, their difference, the monotonicity of the ratio over the whole sweep, and the exact equality of the diagonal row. The last of those is the one most likely to break silently: a change that made this operator’s diagonal non-constant would leave every other number in the essay standing and quietly turn an identity into an approximation.

The derivative that is only half accurate #

The commonest matrix-free operator in nonlinear work is a difference quotient. A Newton method needs J(x)v and can get it without ever forming J:

J(x)v  ≈  ( F(x + εv) − F(x) ) / ε

one extra residual evaluation. Its error has two parts that pull opposite ways: truncation, which is proportional to ε, and cancellation in the subtraction, which is proportional to u/ε. So there is a floor, at ε around √u, and no choice of ε gets below it.

Measured on the Bratu problem — chosen because its Jacobian is a formula, so the quotient is compared against a derivative that is a theorem rather than against a better quotient:

forward difference:  floor 1.3·10⁻¹⁰ at ε = 10⁻⁶,  truncation slope 1.01,  cancellation slope −1.00
central difference:  floor 1.1·10⁻¹²  at ε = 10⁻⁴,  truncation slope 2.00

Both slopes are the theory, fitted rather than asserted. And the floor is six orders of magnitude above the analytic derivative, which is exact.

And it costs nothing #

Here is the result this essay was written for, and it is the counterweight the whole phase needs.

Newton on the Bratu problem, with the analytic Jacobian and with the differenced one:

analytic:            8.00   5.199·10⁻²   3.014·10⁻⁶   5.17·10⁻¹³   4.50·10⁻¹³
differenced, 10⁻⁶:   8.00   5.199·10⁻²   3.010·10⁻⁶   4.59·10⁻¹³   3.87·10⁻¹³
differenced, 10⁻⁸:   8.00   5.199·10⁻²   3.014·10⁻⁶   5.09·10⁻¹³   5.46·10⁻¹³

The same sequence to three digits, and the same final residual. A Jacobian wrong in the tenth decimal place produces a Newton method indistinguishable from one with an exact derivative.

The reason is worth stating exactly, because it is a general principle and not a fact about this problem. The derivative appears only in the step, and the residual is evaluated exactly. An error in J moves where the next iterate goes; it does not move where the fixed point is, because the fixed point is defined by F(x) = 0 and F is computed to full precision. A wrong derivative costs iterations and buys nothing in error.

At the crudest step size on the slider — ε = 10⁻², where the Jacobian is wrong in the fourth digit — the entire cost is two extra iterations. From 10⁻⁴ down it is none.

The fully matrix-free version #

The differenced Jacobian above was assembled column by column, so that the only difference between the two runs was the derivative. A real code does not assemble it. It hands the quotient to GMRES as a closure and never builds anything:

op.apply = (v) => (F(x + εv) − F(x)) / ε

Jacobian-free Newton–Krylov. On the same problem it lands 1.4·10⁻¹⁵ from the analytic Newton answer, in 448 products, with no array of numbers existing anywhere in the computation.

That number — 448 products for eight Newton steps — is where the cost has gone. It is the preconditioning problem from the top of this essay, arriving in its usual form: the inner Krylov solve needs more iterations than it would with a preconditioner nobody can build.

Choosing ε, which is a scaling question #

The difference quotient has one parameter and choosing it badly is the commonest defect in matrix-free code, so it is worth writing the rule down.

The error is roughly ε‖J″‖/2 + 2u‖F‖/(ε‖v‖), minimised at

ε  ≈  2√( u ‖F‖ / (‖v‖ ‖J″‖) )

and the practical version of that — the one in every JFNK implementation — replaces the unknown second derivative by a scale built from the current iterate:

ε  =  √u · (1 + ‖x‖) / ‖v‖

The (1 + ‖x‖) is the important part and it is what a naive implementation leaves out. A fixed ε = 10⁻⁷ is a relative perturbation of 10⁻⁷ only when x is of order one; on a problem whose solution is of order 10⁶ it is a relative perturbation of 10⁻¹³ and the quotient is pure cancellation, and on one of order 10⁻⁶ it is a relative perturbation of 10⁻¹ and the quotient is pure truncation.

That is the scaling argument the previous phase made about condition numbers, arriving in a place nobody expects it: a constant with units in it is a constant about the problem’s description, and a differencing step has the units of x.

What a matrix-free code can still precondition with #

Not nothing, and the list is worth having because it is the practical content of the whole essay.

Physics-based preconditioners. A simplified operator that can be assembled — a lower-order discretisation, a scalar diffusion in place of a tensor one, a linearisation about a simpler state. This is what most production matrix-free codes actually do.

Multigrid. A hierarchy needs a coarse operator, and the geometric version builds it from the grid rather than from the matrix. That is why geometric multigrid is the preconditioner of choice for matrix-free discretisations, and why the algebraic version — which builds the hierarchy from the matrix entries — is exactly the one that is unavailable.

Sketching. A randomised preconditioner is built from products with random vectors, which is precisely the available interface, and it is what changing the condition number on purpose does without ever reading an entry. It is one of the very few genuinely general options.

The diagonal, where it varies. Free here only because this operator’s diagonal is constant; on a variable-coefficient problem it is worth having and costs n probes.

What this phase found, in one place #

Eleven essays have replaced an object by an action: a determinant by an accumulated logarithm, an inverse by a pair of solves, an update by a refactorisation, an eigendecomposition by a polynomial in A, an n×n exponential by an m×m one, an SVD by an iteration, an n²×n² system by a back-substitution. In every case the substitute was better rather than merely cheaper.

This essay is the exception that gives the rule its shape, and it does so twice.

The substitution that costs nothing. A derivative six orders of magnitude worse than the exact one produces the identical answer, because the object it feeds into is not the object the answer is defined by.

The substitution that costs something structural. Giving up the entries costs no accuracy at all and costs a factor of 2.7 in iterations — and what was lost was not precision but every algorithm that needed to look. The tax is on the method library, not on the arithmetic.

So the phase’s sentence needs its second half. The object a question names is usually not the object worth computing; and where the substitution costs something, what it costs is access rather than accuracy.

Where the operators come from, and what each can supply #

Matrix-free is one word covering four situations, and they differ in what they can be asked for besides a product. It is worth separating them, because the preconditioning question above has a different answer in each.

A discretised differential operator. A finite-element or finite-difference code applies the operator element by element or stencil by stencil. It can usually supply the diagonal cheaply — every element contributes to it — and it can supply a coarser version of itself, because the discretisation is parameterised by a mesh. That is why geometric multigrid is the preconditioner of choice here: the hierarchy comes from the modelling rather than from the matrix.

A transform-based operator. A spectral method applies A by transforming, multiplying by a diagonal, and transforming back. Its diagonal in the transform basis is explicit and its diagonal in the physical basis is not. Preconditioning is done in the transform basis, where the operator is diagonal and the problem is trivial — which is the whole reason to use one.

A differenced Jacobian. The case this essay measured. It can supply nothing but products, not even its own diagonal without n residual evaluations, and its products are accurate to about √u rather than u. It is the least informative operator of the four and the most common in nonlinear work.

A kernel or integral operator. A dense matrix that is never formed because it is n² entries of a smooth function. It can supply any entry on demand, which makes it the most informative of the four: hierarchical and low-rank methods work precisely because a block far from the diagonal is numerically low rank and can be sampled.

The list is worth having because it inverts the usual framing. Matrix-free is not a single loss of information; it is a spectrum, and the question for a given code is which of the entries it could supply if asked, and at what price. A differenced Jacobian’s diagonal costs n residual evaluations; a kernel’s costs nothing. The essay’s factor of 2.7 is the price of assuming the worst case when the answer might have been cheap.

That reframing suggests a question worth asking of any matrix-free code before accepting the tax: what would it cost to build the sparsity pattern alone, without the values? For a finite-element operator that is the element connectivity, which the code already has; for a stencil it is the stencil. And a pattern without values is enough for a fill-reducing ordering, for a graph-based coarsening, and for the structure of an incomplete factorisation — leaving only the values to be filled in, at one probe per colour of a graph colouring rather than one per column. Probing a sparse operator for its entries costs the number of colours, not n, and on a stencil that is a small constant. The factor of 2.7 above is the price of never asking.

What is worth carrying #

A Krylov method needs one thing from a matrix and it is not the entries. Everything built on products survives a matrix that does not exist; everything built on comparisons does not.

What is lost is the method library, and it costs iterations rather than digits. A factor of 2.7 on the model problem, growing with the size, because the preconditioner that works is defined by a pattern.

A differencing step has the units of x, so a fixed ε is a bug that only appears on problems scaled differently from the one it was tuned on. Scale it by (1 + ‖x‖)/‖v‖ and it is a relative perturbation everywhere.

And a derivative wrong in the tenth digit costs nothing, because it appears only in the step while the residual is evaluated exactly. The general form: a method whose answer is defined by a quantity it computes exactly can afford to be very wrong about everything else — and it is worth asking, of any approximation, which of those two roles it is playing.