A hierarchy with no grid behind it

The previous essay establishes that the construction does what it says: given a matrix and nothing else it reproduces the geometric operators on the model problem and discovers semi-coarsening on the anisotropic one. That is the interesting half of the claim and it is the half that gets written down.

This essay is the other half. What does the hierarchy cost, on a problem there is no geometric alternative for, and where does it stop?

The test matrix is a graph Laplacian: 400 vertices arranged as a ring with about 1,200 random chords, weighted, with a small diagonal shift to make it positive definite. Nothing is anywhere. The couplings have no direction, the degrees run from three to about a dozen, and there is no coarse grid to write down because there is no grid. Every component in the geometric essays is inapplicable by construction — not worse, inapplicable.

The algebraic V-cycle converges on it at 0.199 a cycle, which is a rate a Krylov method would be pleased with — and which is not the best this hierarchy can do, for a reason two sections below. The rates the multigrid field is used to quoting come from a rate that does not notice the size, where the operator is a discretisation and every parameter has been calibrated on one. And it costs a great deal more than that number suggests.

Two complexities, and the one that gets quoted #

A hierarchy is measured by two ratios and they are routinely treated as one.

Grid complexity is the total number of unknowns over all levels, divided by the fine grid’s. It says how much extra storage for solution vectors the hierarchy needs, and how much of the smoothing work happens below the top level.

Operator complexity is the total number of stored matrix entries over all levels, divided by the fine matrix’s. It says how much extra storage for operators is needed, and — because a smoothing sweep costs one multiply-add per stored entry — how much arithmetic a cycle actually performs.

On the geometric hierarchies of this phase the two are close: 1.296 and 1.483 for full coarsening on the 31×31 grid. Nothing about the sparsity changes down that hierarchy, so counting unknowns and counting entries give nearly the same answer, and it is easy to form the habit of quoting either.

On this graph Laplacian they are 3.05 and 17.7.

The grid complexity is unremarkable. A reader given only that number would conclude the hierarchy costs three times a single-level method, which is what a geometric hierarchy in three dimensions costs and nobody objects to. The operator complexity says a V-cycle performs nearly eighteen times the arithmetic of one fine-grid sweep, which is a completely different method.

Where the entries come from #

They come from the Galerkin product, and the mechanism is the one the sparsity field spent three essays on: eliminating a variable couples everything it touched to everything else it touched.

PᵀAP does exactly that, in one step, for every dropped point at once. A coarse point that interpolates to several fine points inherits couplings to every coarse point those fine points touch, so the coarse operator’s graph is roughly the square of the fine one’s restricted to the coarse set. On a grid that squaring is contained: a five-point stencil becomes a nine-point one and then stops, because the geometry limits how many neighbours a neighbour can have. On a graph with random chords there is no such limit, and the squaring compounds level after level.

The per-level counts say it plainly. The fine matrix has 400 unknowns and 2,784 stored entries — seven a row. Level 2 has 185 unknowns and 10,989 entries: 59 a row. Level 3 has 135 unknowns and 14,969 entries, which is 111 a row on a matrix with 135 columns, so that level is 82% dense. And level 6, with 41 unknowns, stores all 1,681 of its entries: it is a dense matrix.

A hierarchy that removes a third of the unknowns per level and squares the couplings has not made the problem smaller in the way a count of unknowns suggests. The two complexities disagreeing by a factor of six is the measurement that says so, and only one of them would have said it.

The rate degrades with the density too #

The cost and the rate move together, which is worth measuring because a reader might reasonably expect a denser hierarchy to be a more accurate one.

Take the same construction and vary how many chords the graph has:

Past about one chord a vertex the method gets worse in both currencies at once. The coarse operators are denser and the corrections they produce are less useful, because interpolation from a neighbourhood that is nearly everything is interpolation that says nothing in particular.

The table’s step is too coarse to find its own minimum, though, and filling in the half-steps moves it.

So the optimum is at 1.5 rather than 1.0, and it is a factor of two below the row the table calls best. Between 1.5 and 2.0 the factor then jumps from 0.030 to 0.199 — nearly sevenfold across half a chord — so the cliff the section is about is both steeper and later than five rows can show.

Past the cliff the damage stops accelerating. From 2.5 to 3.5 chords the factor costs 0.233; from 3.5 to 4.0 it costs 0.067 across half a chord, which is 0.134 a chord — a little over half the earlier rate. The curve is bending towards a ceiling rather than running away, so a hierarchy built well past its optimum is wasteful rather than useless.

And the complexity that gets quoted is the one that does not move. Across the whole sweep the grid complexity runs 2.56, 2.78, 3.06, 3.11, 3.01, 2.99 — a spread of 21% — while the convergence factor goes from 0.030 to 0.671, a factor of twenty-two. A reader handed the grid complexity would see a hierarchy whose shape barely changed; the method meanwhile went from converging in a handful of cycles to barely converging at all. That is the same point the earlier section made about which number gets reported, arriving with a measurement behind it: the quoted complexity is not merely incomplete, it is uncorrelated with the thing a caller cares about over the range where the parameter actually does damage.

That reads as a statement about the method’s applicability, and it is a statement about a default.

The same sweep with the threshold chosen for the matrix #

Every number in that table was taken at the conventional θ = 0.25, and the section below on the threshold shows that 0.25 is several times worse than the best available on a graph like this. Doing the scan at every chord density instead:

The degradation is the threshold. At three chords a vertex the default gives 0.492 and θ = 0.1 gives 0.0161 — thirty times better, and better than the sparsest graph in the sweep manages at any threshold. The best rate anywhere in the sweep belongs to the second-densest graph, not to a sparse one.

And the retuning improves the storage as well: 17.7 → 12.3, 19.6 → 15.4, 18.2 → 16.5. So worse in both currencies at once is exactly inverted — at a threshold chosen for the matrix, a denser graph is better in both.

The direction of the retuning is the part worth carrying. The best θ falls as the graph gets denser, from 0.20 to 0.05: what a dense graph wants is more of its couplings called strong, not fewer. That is not the intuition — a graph with many chords sounds like one where a strength test should be more selective — and it is what makes 0.25 a poor default here rather than a mildly suboptimal one.

So the honest shape of classical algebraic multigrid’s applicability is narrower than the paragraph above claims. There is no boundary at one chord a vertex. What there is is a threshold calibrated on discretised operators, applied unchanged to a matrix that is not one, and the cost of that is a factor of thirty in the rate and five in the storage — which is an actionable finding where an applicability boundary is only a warning.

There is a caution attached to it, and it is the one choosing without knowing exists for. The best θ in that table was found by solving the problem eight times, which is not a procedure anybody runs before solving the problem once. What the measurement licenses is not tune θ — it is do not read a bad rate at θ = 0.25 as a fact about the matrix, because on this family it is a fact about the default. A code that sees a poor factor on an unstructured matrix should try two or three thresholds before concluding anything, and the two or three that are worth trying are below 0.25 rather than above it. The same shape of advice appears in the sparsity field, where the order decides the memory finds four orderings within a factor of 1.7 of each other and one default that is not among them.

The other method in this field pays the same threshold, which is why the finding travels. Smoothed aggregation reads the identical strength graph, and its own failure at 45° is diagnosed there as a property of the matrix rather than of θ — a diagnosis that would be worth re-examining in the same way.

What the second pass does here #

The mechanism essay reports that the second coarsening pass — the one that promotes an F-point when two F-points depend on each other with no shared C-point — does nothing at all on the one-dimensional Laplacian: zero promotions.

On this graph it supplies 56% of the coarse set. Of the 277 points kept at the top level, 155 were promoted by the second pass rather than chosen by the first.

That is the clearest single indicator of how different an unstructured matrix is. On a grid, the first pass’s greedy choice already produces a set from which everything can be interpolated, because the neighbourhoods overlap regularly. On a graph they do not, and more than half of what is kept is kept because the interpolation would otherwise have had nothing to work from. It is also why the coarsening ratio is 0.69 rather than the grid’s 0.5 — the method keeps what it must, and here it must keep a great deal.

Where the geometric method applies, the algebraic one is dearer #

The graph is the case with no alternative. The isotropic model problem is the case with one, and it is the fair comparison: the same 31×31 operator, the same Galerkin products, and the only difference the choice of which points to keep and how to weight them.

The stored entries per row are nearly the same at every level — 8.22 for the geometric first coarse operator against 8.48 for the algebraic one, both of them the nine-point stencil the Galerkin product produces from a five-point one. What differs is how fast the unknowns fall away. Geometric full coarsening quarters them; the algebraic coarsening, choosing from the entries, keeps 481 of 961 — just over half — and so needs a level whose operator is four times the size of the geometric one’s.

That is not a defect either. The algebraic construction is choosing a maximal set of points no two of which are strongly connected, which on the five-point stencil is a chequerboard, and a chequerboard is half the points rather than a quarter. Geometric coarsening takes a quarter because somebody knew the grid was two-dimensional and could halve both directions at once. The algebraic method has one relation and takes what it implies.

Both converge: 0.0656 for the algebraic hierarchy against 0.2016 for the geometric one — the algebraic method being the faster of the two per cycle, at 1.47 times the arithmetic, which is close enough to a wash that neither is worth preferring on this problem. The four extra unknowns per coarse point are the price of not being told what the grid was.

The setup is a third cost, and neither complexity measures it #

Both ratios above describe a cycle. Before any cycle runs, the hierarchy has to be built: the strength sets computed for every level, the two coarsening passes run, the interpolation assembled, and the triple product PᵀAP formed at each level.

The triple product is the expensive part, and it is expensive for the same reason the coarse operators are dense: forming it touches every path of length two through the fine matrix restricted to the coarse set. A geometric hierarchy avoids all of this — its transfer operators are three numbers known in advance and its coarse operators can be rediscretised — so the setup is free there and is not free here.

What that means practically is that the two complexities understate the algebraic method’s cost on a matrix solved once, and overstate it on a matrix solved many times with the same structure, where the setup is amortised over every solve. A time-stepping code that rebuilds the same operator every step is in the first case; a Newton iteration reusing a hierarchy across several steps is in the second. This phase measures neither — every number here is a convergence factor or a stored entry, and the setup cost is a real omission rather than a hidden one.

The threshold is a real parameter, and 0.25 is not the right value #

On the anisotropic grid operator the strength threshold has a meaning and the answer is insensitive to it either side of ε. On an unstructured graph there is no ε and no such interval, and the scan is not reassuring:

The factor moves by a factor of 43 across the range, the ordering is not monotone, and θ = 0.25 — the value every introduction gives — is ten times worse than θ = 0.20 on this matrix. That is not an argument for a different constant; θ = 0.2 is not better in general either. It is the measurement that says the constant is doing less work than its ubiquity suggests, on problems away from the grids it was calibrated on.

There is a trade underneath the noise, and it runs the useful way. Among the thresholds that converge as fast as the best, θ = 0.05 stores 8.05 where the conventional value stores 13.44 — 40% less arithmetic per cycle for the same rate. Nothing in the method finds that; a scan does.

Where it stops: forty-five degrees #

The rotated anisotropic operator is where this phase’s iterative half ends, and the ending is not a success.

At 45° with ε = 10⁻³ on a 31×31 grid:

The algebraic method is the worse of the two families here — worse than the geometric method it was brought in to replace. That is not the shape of result a phase goes looking for and it is the one the measurement gives.

The reason is locatable and it is not a defect of the algebraic construction. At 45° the strong direction of the equation runs along the anti-diagonal. The standard nine-point discretisation of that operator has couplings of −0.5005 along each grid axis and −0.2498 along the anti-diagonal, so the matrix’s largest entries are not in the equation’s strong direction — the discretisation has smeared the anisotropy across the axes. It also has two positive off-diagonal entries, which the strength measure correctly ignores and which put the matrix outside every M-matrix theory the field has.

An algebraic method sees the matrix. It cannot recover an anisotropy the discretisation did not encode, and no threshold rescues that: the strength measure is reading real entries and reporting them correctly.

So the honest statement is that the difficulty at 45° is partly a discretisation problem rather than a solver problem, and that this phase measured three solvers against it and none of them worked. A better discretisation of the rotated operator, or an interpolation built by minimising energy rather than by distributing row weight, are the two directions out, and neither is measured here.

What is asserted, and what refuses #

It converges on a matrix with no geometry: 0.199 a cycle, at least five levels, on 400 unknowns and 1,192 edges.

The two complexities disagree, asserted as a ratio — operator complexity more than five times grid complexity — rather than as two numbers a reader has to compare.

Some level is more than half dense, asserted by scanning the hierarchy rather than by naming level 6, so that a change in the construction cannot make the assertion vacuous.

The second pass supplies more than 30% of the coarse set, which is the measurement that distinguishes this matrix from a grid.

The factor degrades with density, asserted in both directions: below 0.15 at half a chord a vertex, above 0.5 at four.

The conventional threshold is not the best, asserted as conventional > 5 × best — a check that would fail if some future change made θ = 0.25 respectable, which is the point of writing it that way.

And a hierarchy that never coarsened is refused. At θ above one nothing is strong, no point is dropped, and the “V-cycle” is a dense direct solve on the fine matrix reporting a convergence factor of zero. The solve throws rather than returning it.

What is left #

Smoothed aggregation, which builds coarse variables by grouping fine ones and applies a smoothing step to the resulting interpolation. Its operator complexities are typically far lower than classical AMG’s on exactly the matrices this essay struggles with, and it is not measured here.

Energy-minimising interpolation, which is the standard answer to the rotated case and to problems where distributing a row’s weight is the wrong prescription.

Aggressive coarsening, which drops more points per level to hold the complexity down and needs a longer-range interpolation to compensate. It is the direct response to the numbers in this essay and it is a body of work rather than a variation.

What is settled is the measurement the phase set out to make. An algebraic hierarchy is not a free replacement for a geometric one. It is applicable where nothing else is, it reproduces the geometric answer where the geometric answer exists, it costs between one and a half and twelve times as much arithmetic per cycle depending on the matrix, and on the one problem all three methods were run against it came third.

The hierarchy a matrix-free code cannot build #

This page’s hierarchy comes out of the matrix entries and nothing else, which is exactly what makes it unavailable where the matrix is a subroutine. The geometric version, which reads the grid instead, is the one matrix-free codes use.