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.

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.

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 depth phase’s line smoother and semi-coarsening return 0.967 on, and smoothed aggregation returns 0.789 on.

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