The formula that was already optimal

This site has built a coarse space three ways.

The classical algebraic hierarchy derives each interpolation weight from one row of the matrix, by assuming the error is smooth along strong connections and solving that row for the fine value. Smoothed aggregation partitions the points, builds a piecewise-constant prolongator, and applies one Jacobi sweep to its columns.

The third way asks a question instead of following a recipe: which P, over a given sparsity pattern, minimises the energy of its own columns?

minimise Σⱼ PⱼᵀAPⱼ, optionally subject to P reproducing a near-null-space vector exactly

Written out that is a quadratic in the nonzeros of P with linear equality constraints — one saddle-point system, solved once. It is the deferral the depth phase wrote down and the synthesis phase left — one saddle-point system, solved once.

That is drawn at the strongest anisotropy the figure offers. Turning it down says which parts of the picture are about the rotation and which are about the anisotropy.

The gap appears all at once and then stops growing, which is not what a reader would guess from the first frame.

The energy gap converges to about 1.203 rather than growing. 1.0000, 1.2230, 1.2047, 1.2029 — it appears between ε = 1 and ε = 0.1, overshoots slightly, and settles. So the classical formula is twenty per cent above the least energy available at 45°, whatever the anisotropy, and there is no regime in which it becomes arbitrarily bad. That is a much more useful statement than “the two part company”: it says how far apart they get and that the answer is bounded.

At 0° the ratio is one at every anisotropy, to four decimal places at three of the four stops. The classical formula is exactly optimal along the grid axis however extreme the operator, which is what makes the 45° reading a diagnostic rather than a general defect — a code seeing a ratio of one has learnt that its strong direction is aligned, and that is knowledge it had no other way to get.

And the ordering of the three methods inverts. At ε = 1 the constrained variant is the worst of the three, 0.0848 against 0.0609 for both others; by ε = 0.001 it is the best, 0.2131 against 0.2843 and 0.3464. The constraint costs something when there is no anisotropy to exploit and pays when there is, so it is not a strict improvement and a code that applies it unconditionally is paying on the isotropic problems to win on the rotated ones.

One more reading, because the two quantities are not the same size: a 20% energy gap buys a 38% better rate — 0.3464 down to 0.2131 at the finest ε. Energy is what the minimisation optimises and the rate is what a caller waits for, and the second responds nearly twice as strongly as the first.

The first answer is that there is no third way #

On the one-dimensional Laplacian, with the coarsening the classical method chooses, the energy minimiser and the classical formula produce the same matrix. Not to a tolerance — the worst difference over every row and every column is zero, and both reach an energy of 15.000000 and a two-level convergence factor of 0.06126.

Two derivations with nothing in common. One is local, algebraic and row-by-row: take row i of A, assume the error at the fine point is a weighted average of its strong neighbours, solve. The other is global and variational: 30 unknowns, a quadratic form, a linear solve of the whole thing at once.

That is a two routes to a number result of an unusual kind, because the two routes are not two arithmetics for one formula — they are two derivations that were never supposed to meet.

It survives into two dimensions on the operators the row-wise assumption is right about:

Energies agreeing to six digits and convergence factors to three.

The linear term, which the first version dropped #

The objective has a piece that is easy to leave out, and leaving it out produces a coherent, wrong, and interesting-looking answer — which is worth recording because it took a check to notice.

A coarse point interpolates itself, so its row of P is the identity and is not a variable. That means every column of P already contains a one before any variable is chosen, and the energy has cross terms between that fixed entry and every variable in the same column: the objective is a quadratic plus a linear term, not a pure quadratic.

Dropped, the free minimum of a positive-definite quadratic with no linear term is P = 0. Every fine row comes back empty, the method looks useless, and the constrained version looks like the only one worth having — which is exactly the story one expects to find, and is why it survived a first reading — the same shape as a plausible number hiding a defect.

With the linear term in place the free minimum is the classical interpolation, and the whole essay changes shape.

Where they part company #

Rotate the anisotropic operator so its strong direction runs at 45° to the grid — the problem the line smoother and semi-coarsening return 0.967 on, and smoothed aggregation returns 0.789 on.

And in between, the minimiser is worse #

Nought degrees and forty-five are the ends of an axis, and the middle of it is not a line between them. Sweeping the rotation on the same operator:

The energy column behaves exactly as it must: the minimiser is at or below the classical formula at every angle, because that is what it minimises.

The rate column does not. From 5° to 22.5° the classical formula converges faster — by a factor of 1.6 at ten degrees, where the minimiser’s energy is 1.3 per cent lower and its convergence factor is 0.1420 against 0.0889. Four of the seven angles measured go that way.

And the two gaps do not track each other even in sign. The largest energy gap is at 30°, where the minimiser has 27 per cent less energy and the rates are within three per cent of each other. The largest rate gap is at 45°, where the energy gap is smaller.

Which says what the objective is worth #

Two cautions on reading that table, both of which narrow it rather than dissolve it.

The rates are two-level convergence factors on one grid size, and a two-level factor is not what a V-cycle delivers — the recursion compounds it, and a construction that is 1.6 times worse per level is worse than that over a hierarchy. So the reversal is if anything understated by these numbers rather than an artefact of them.

And the sweep varies the angle at one anisotropy, ε = 10⁻³. The angle and the strength are two knobs and only one is moving here; a milder anisotropy would move the whole reversal, and where it moves it to is not measured. What the sweep establishes is that the reversal exists and is not confined to a single angle, which is what the two endpoints could not say.

That is not a defect in the construction. It is the objective being a proxy, and the measurement is how good a proxy it is.

Column energy is chosen because minimising it is a solvable problem — one saddle-point system, solved once — and because a low-energy interpolation is, loosely, one that represents the operator’s smooth modes well. The quantity that actually decides the method is the two-level convergence factor, which is a spectral radius of a product of projections and is not a quadratic in anything. So the construction optimises the tractable quantity and hopes for the other, which is the ordinary situation and is worth being explicit about rather than assumed.

What the sweep adds is that the hope fails over a third of the axis, and fails in the region where the classical assumption is starting to break but has not broken. At 0° the row-wise formula is exactly right and there is nothing to gain. At 45° it is wrong and the minimiser gains a factor of 1.5. Between them the classical formula is approximately right and the minimiser’s extra freedom is spent buying energy the convergence does not want.

That gives the practical rule a shape the two endpoints alone do not. It is not use the minimiser where the operator is not aligned. It is: the minimiser pays where the row-wise assumption fails outright, and costs where it merely bends. Which is harder to detect than either endpoint suggests, because the quantity that says which regime a problem is in is the convergence factor — the thing being predicted.

The honest reading of the whole essay is then two findings rather than one, and both are about the hierarchy the matrix chooses rather than about a grid. On the operators the row-wise assumption is right about, the two constructions are the same matrix, which is the surprise. On the operators it is wrong about, the minimiser is better. And on the operators in between — neither, and the objective points the wrong way.

| minimiser, constrained to reproduce a constant | 103.27 | 0.2600 |

The classical interpolation is 22% above the minimum energy over its own sparsity pattern, and converges at 0.527 where the minimiser converges at 0.345.

So the agreement in the first half of this essay is not a theorem. It is a statement about the operators the row-wise assumption is correct for, and the assumption is exactly the one the rotated stencil breaks: the error is smooth along strong connections, and the rotated stencil’s axis couplings are −0.5005 against a diagonal coupling of −0.2498, with two of the four diagonal entries positive — which no strength measure on this site reads correctly.

The energy gap is a diagnostic #

That is the practically useful part, and it is worth separating from the comparison.

Computing the energy of a given P costs a triple product and no reference solution. Computing the minimum over the same pattern costs one saddle-point solve. The ratio of the two is then a number that says whether the interpolation in hand is the best one available on the pattern it occupies — and it needs no exact answer, no known spectrum, and no second method to compare against.

A ratio of one says the classical formula has nothing left to give on this pattern; the way to improve is to widen the pattern. A ratio of 1.22 says the formula’s derivation has stopped applying and a different P on the same pattern would do better.

Every other diagnosis of this failure on this site has needed something extra: the depth phase read it off the stencil coefficients by hand, the synthesis phase predicted it from a strength measure and then checked it against a solve. This one is a quotient of two numbers computed from A and P.

The pattern is the other half, and it is the bigger half #

Every number above is a minimum over a fixed sparsity pattern — the one the classical formula occupies, which is each fine point’s coarse neighbours. That pattern is a decision too, and it is usually made without comment.

Widen it: allow each fine row to reach the coarse points two steps away in the matrix graph. On an 11×11 grid that takes P from 165 nonzeros to 365.

On the aligned operator, doubling the storage buys nothing — the same energy to six digits and the same convergence factor. The optimum over the narrow pattern is already the optimum over the wide one, which is a sharper statement than the classical formula is good enough: there is no interpolation on twice the pattern that does better.

On the rotated operator it buys a factor of 1.8 in the convergence factor.

And it buys that for a 2% fall in energy. That pair is the honest limit of the diagnostic proposed above: energy and convergence move together in sign and not in size, so a small energy gap can hide a large difference in what the method does. The ratio says something is wrong here; it does not say how much.

What it would cost to do this for real #

The measurements above solve a dense saddle-point system of size (nonzeros of P) + (fine rows) — 702 unknowns for a 225-point grid, and cubic in that. Nobody does this.

What the practical methods do instead is take a few steps of a Krylov method on the constrained problem, starting from the classical interpolation, and stop long before convergence. That is worth noticing in the light of the first finding: if the classical interpolation is the unconstrained minimiser, then starting there and taking a few constrained steps is a method for imposing the constraint rather than for minimising anything — the minimisation is already done before the first step.

So the cost of the exact answer is what makes the exact answer worth computing once, on a small problem, as a measurement rather than as a method. Which is the position this whole site takes about its exact ground truths: exact.js solves in BigInt rationals nobody would use in production, so that the float answer has something to be wrong against.

The constraint is the modelling #

The near-null-space constraint — every row of P sums to one, so a constant is interpolated exactly — is what energy-minimising interpolation is usually described as being for. Measured on its own, it is a decision with a price in both directions.

On the two problems the unconstrained minimiser handles well, forcing the constant makes the method up to five times worse. On the one it handles badly, the same constraint is worth a third.

The cost has a location. These are Dirichlet problems, so the constant is not in the operator’s near-null space: the rows next to the boundary do not sum to zero, and forcing P to reproduce a constant there forces the wrong vector. Fifty-six rows of 225 in two dimensions, and that is enough to cost a factor of five.

And the benefit has a location too. On the rotated operator the near-null space is precisely what the classical formula has lost, so an assumption that puts it back — even the wrong one, even at the boundary — is worth more than it costs.

What this says about the near-null space #

The reading the measurements support is sharper than supply a good near-null-space vector.

The near-null-space constraint is not an approximation to be improved. It is a statement about the problem, it is either right or wrong, and where it is wrong it costs a factor of five while looking exactly as principled as where it is right.

And the classical formula already encodes one. Its row-wise derivation — the error is smooth along strong connections — is a near-null-space assumption written in a different notation, and where that assumption holds the formula is optimal without any constraint being imposed. The two methods are not a heuristic and an optimisation; they are two ways of stating the same modelling, and they differ exactly where the statements differ.

What was actually solved #

The minimisation is a saddle-point system and it is worth saying exactly which one, because the answer depends on choices that are easy to leave implicit.

The variables are the nonzeros of P’s fine rows. A coarse point interpolates itself, so its row is the identity and is fixed — not a variable, not a constraint, simply given. The objective is Σⱼ PⱼᵀAPⱼ with P assembled from the fixed rows and the variables, which makes it a quadratic plus a linear term in the variables. The constraints, when present, are one equality per fine row.

Two of those choices are load-bearing.

Fixing the coarse rows is what makes the problem well posed. Left free, they minimise too, and the free minimum sets them to zero — a P with no identity in it, a singular coarse operator, and a two-level method that cannot be run at all.

And the linear term is what makes the free minimum interesting. Without it the objective is a pure positive-definite quadratic whose minimum is P = 0, which is a coherent answer to a different question. The first version of this measurement asked that question, got that answer, and produced an essay in which the coarse problem was never named as which energy minimisation only worked when constrained. Both the code and the conclusion were self-consistent and wrong.

The fourth method at 45° #

This phase has now put four constructions on the same rotated operator, and the sequence is worth laying out because the last entry breaks it.

Convergence factors, from three phases of this site and this essay.

The first four have the same shape: excellent aligned, hopeless rotated. The fifth is the only one whose two columns are close together — and it gets there by being worse on the aligned problem, not by being better on the rotated one in isolation.

That is not a disappointing result. It is the trade stated plainly: the aligned methods are aligned methods, they encode a direction, and their aligned numbers are what the encoding buys. A method that encodes less does worse where the encoding was right and better where it was wrong, and the whole of what this essay measures is where the boundary between those two regions is.

What no method in the table does is be best in both columns. The nearest thing to a general recommendation the measurements support is that the diagnostic generalises even though the methods do not — the energy ratio is computable for any of them, on any operator, and it is one exactly when the method in hand has nothing left to give on the pattern it occupies.

What the drag does #

The slider is the anisotropy, from a ratio of one to a thousand.

At a ratio of one the operator has no strong direction to lose, so rotating it changes nothing: the energy ratio is 1.0000 at every angle, and every interpolation converges at about 0.06. That is the negative control and it is a position of the slider rather than a sentence — the assertion is written in two halves so that the figure at that setting is required to show no gap.

As the anisotropy sharpens, the rotation has more to spoil, and the gap at 45° grows to 22%. What does not change at any position is the value at zero: the classical formula is the minimiser on an aligned operator whatever the anisotropy ratio, which is the half of this essay that is not about failure.

A formula that is optimal and unusable #

A closed form that turns out to be the minimiser of something is the happiest result there is. The unhappy version is a closed form that is exactly correct and is not a computation of what it names.