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.

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. 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 grows with the problem: 1.93 at 64 unknowns, 2.67 at 400. So matrix-free is not a fixed tax; it is a tax that increases with size — which is awkward, because the reason the entries are unavailable is usually that the problem is large.

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. 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.