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

That is the honest shape of classical algebraic multigrid’s applicability. It is built for matrices whose couplings are local in some sense the strength measure can find, and a graph with enough long-range chords stops being such a matrix without ever stopping being sparse.

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.