The coarse grid the matrix chooses

Both repairs for anisotropy work and both are aimed. A y-line smoother and a y-semi-coarsening each take the convergence factor from 0.957 to under 0.11, and each returns 0.966 on the same problem with the axes swapped. The direction was written into the code by whoever typed the string; nothing read it off the operator.

Everything the geometric method needed was a question about couplings. Which points can be dropped? The ones whose values are determined by their neighbours. How should a dropped point be recovered? From the neighbours that determine it. What is the coarse problem? Whatever the fine operator looks like from the coarse points. And the couplings are in the matrix.

Algebraic multigrid answers all three from the entries of A and nothing else. This essay is the mechanism and the checks that it is doing what it claims; the next one is what it costs and where it stops.

The claim, and how it is kept honest #

The claim is that the hierarchy comes out of the matrix. A claim of that shape is easy to make and easy to violate by accident — one import of a grid helper, one length that happens to be the square root of the dimension, and the code is quietly geometric again.

So the arrangement here is structural. The algebraic construction lives in a file that imports from one place: a sparse matrix type that knows about rows and columns and nothing about space. It does not import from the geometric file, and the geometric file imports from it for the comparisons, so the dependency runs one way. The claim is checkable by reading an import list, which is worth more than a paragraph promising the same thing.

Where a grid appears in this essay it is there so a reader can see what the construction arrived at. The construction never had it.

Question one: which couplings matter #

Point i strongly depends on point j when

−a_ij  ≥  θ · max over k ≠ i of (−a_ik)

with θ a threshold, conventionally 0.25. Three things about that definition earn their place.

The minus sign is not decoration. For the operators this subject produces — an M-matrix with a zero row sum, which is what a discretised diffusion problem gives — a large negative off-diagonal means the two values are pulled together, and a positive one means they are pushed apart. Treating a positive coupling as strong produces interpolation weights of the wrong sign. The rotated operator from the last essay has two positive off-diagonals in every interior row, so this is not hypothetical.

It is a ratio within a row, not an absolute size. A row whose couplings are all 10⁻⁶ has strong connections; a row with one coupling of 1 and one of 10⁻³ has one. That is what makes the measure scale-free and what makes it find anisotropy.

It need not be symmetric. i can strongly depend on j while j does not strongly depend on i, whenever the two rows have different largest couplings — which is why the coarsening below keeps both directions and distinguishes “depends on” from “is depended on”.

On the operators in this essay it happens to be symmetric anyway: measured on the anisotropic operator at ε = 0.3, every strong pair agrees, so the fraction is exactly 1.000. That is a property of these matrices rather than of the definition, and it is worth knowing which. A matrix with rows of genuinely different scales — a discretisation on a graded mesh, a system assembled from elements of different sizes — breaks it immediately, and the code that assumed an undirected graph would then be reading a relation that is not there.

What “smooth” means when there is no grid #

The geometric method’s whole vocabulary rests on a picture. An error is smooth when it varies slowly across the grid; a coarse grid can represent it; a smoother removes the rest. Take the grid away and every word of that is undefined, so the construction needs a replacement, and the replacement is one this site has already spent an essay on.

An algebraically smooth error is one the smoother has stopped removing — and a smoother stops removing an error exactly when its residual is small. So the errors the hierarchy has to handle are those with Ae ≈ 0 and e not small: a small residual and a large error, which is the pair a small residual is not a small error exists to separate.

Write out what that means for one row. If (Ae)ᵢ ≈ 0 then

aᵢᵢ eᵢ  ≈  − Σ over j ≠ i of aᵢⱼ eⱼ

so eᵢ is a weighted average of its neighbours’ values, with the weights being the couplings. Where one coupling dominates the row, the average is nearly that neighbour’s value alone: along a strong coupling, an algebraically smooth error varies little. That is the sentence the geometric picture was standing in for, and it says the same thing without mentioning space.

Everything else follows from it. A point whose row is dominated by a few strong couplings can be recovered from those neighbours, so it can be dropped; the weights to recover it with are the couplings themselves, normalised; and “strong” has to be measured as a ratio within the row, because the statement is about which term dominates that particular average.

It also says what the method cannot do. The derivation assumes the residual is small in a sense that makes the row equation a good approximation, which is the smoother’s job to arrange, and a smoother that leaves a rough error — one with a large residual — leaves an error this reasoning says nothing about. The whole construction inherits the assumption the geometric one had, restated in terms the matrix can express.

Question two: which points to keep #

Ruge–Stüben coarsening, in two passes.

The first pass gives every undecided point a measure: how many points strongly depend on it. Take the largest, make it a C-point — kept — and make every undecided point that strongly depends on it an F-point — dropped, to be interpolated. Each new F-point then raises the measure of everything it strongly depends on, because those have become more valuable as interpolation sources. Repeat until nothing is undecided, breaking ties by index so the answer is a function of the matrix and not of the iteration order of a map.

The second pass repairs what the first cannot guarantee. If i and j are both F-points and i strongly depends on j, then j’s value has to reach i somehow, and the only route is a C-point they both depend on. Where there is none, one of the pair is promoted to C.

That second pass is usually described as a refinement and it is not optional: without it the interpolation below meets rows whose strong C-neighbourhood is empty and has to fall back to something arbitrary. How much work it does is a measurement worth reporting, and it varies enormously — zero promotions on the one-dimensional Laplacian, and more than half the coarse set on the unstructured matrices of the next essay.

Question three: how to interpolate #

A C-point is copied. An F-point i is written as a weighted sum over the C-points it strongly depends on, and the weights come from insisting that the whole of row i’s off-diagonal weight is accounted for even though only some of its neighbours are available to carry it:

αᵢ = (sum of all the negative aᵢⱼ) / (sum of the negative aᵢⱼ over kept strong j)
wᵢⱼ = −αᵢ aᵢⱼ / aᵢᵢ

and the same with β for the positive couplings, scaled separately. The two scalings are one line of algebra with a large consequence: when row i sums to zero, the weights sum to one, so interpolating a constant returns that constant.

That is the property the whole method rests on. The error a smoother leaves is locally near-constant along strong couplings — that is what “smooth” means once there is no grid to be smooth on — and an interpolation that could not represent a constant would put it back wrong. It is asserted here on two matrices, on every row whose operator row sums to zero, to 10⁻¹⁰.

The coarse operator is then PᵀAP, the same Galerkin product the geometric hierarchies use. It is symmetric whenever A is, to 10⁻¹⁴ as computed, and it still annihilates a constant away from the boundary — the property the two-dimensional essay found was the one that mattered, once the identity it was thought to rest on turned out to be false.

First check: it rediscovers the geometric answer #

Give it the one-dimensional Laplacian, tridiagonal(−1, 2, −1), on 63 points, and no information about a grid.

Both off-diagonals of every interior row are equally and maximally strong, so the coarsening has to break the tie somehow. What comes out is 31 kept points of 63, no two of them adjacent — alternating, which is what “every other point” means without assuming where the set starts. The second pass promotes nothing. And the interpolation weights it derives are 1 to a kept point’s own copy and ½ to each neighbour, with the largest disagreement against those numbers exactly zero.

Those are the operators the geometric hierarchy was handed on the first page of the multigrid field, where they were justified from a picture of two grids. Here they are consequences of the entries.

Second check: it discovers semi-coarsening #

Give it the anisotropic operator at ε = 10⁻³ on a 31×31 grid.

The x-couplings are 10⁻³ and the y-couplings are 1, so the threshold at θ = 0.25 admits none of the former: zero of the 930 horizontal couplings are strong. Nothing is therefore dropped for depending on its horizontal neighbour, and the coarsening becomes a question about the vertical couplings alone.

What comes out: all 31 grid rows are kept whole or dropped whole, 15 of them are kept, and no two kept rows are adjacent. That is the semi-coarsened grid — the one the previous essay had to be told to build — arrived at without a coordinate, without a mesh spacing, and without anybody naming a direction.

Three numbers rather than one, because keeping half of everything at random would give the same coarsening ratio and a completely different grid.

The threshold is a switch, not a dial #

θ is introduced everywhere as a tuning parameter with a conventional value, which suggests a smooth trade-off and a number worth fiddling with. On this operator it is nothing of the kind.

The x-couplings have size ε and the y-couplings size 1. So θ below ε calls the x-couplings strong and θ above ε does not, and the coarsening switches between full and semi-coarsening at that one point. Either side of it the answer is insensitive to θ over more than a decade; across it the convergence factor moves by an order of magnitude — from about 0.46 to about 0.048 at ε = 10⁻².

So the threshold has a meaning: it is the decision about which couplings count, and on a problem with a scale in it the decision has a right answer with the scale in it. The conventional 0.25 is not a good value or a bad one; it is right for every anisotropy below a quarter and wrong for every one above it.

What it does not have to be told #

The strongest available statement that the direction is not built in is the one the geometric repairs failed: turn the problem through a right angle and measure again.

The algebraic method returns the same number twice, to a thousandth, because the construction it performs on the second matrix is the first construction with the indices relabelled. It is slightly slower than the best geometric repair on the orientation that repair was aimed at, and it does not have an orientation to be aimed at.

That is the whole of what has been bought, and it is worth being precise about what has not. The method is not more accurate, does not converge faster, and is not cheaper — the next essay measures a hierarchy that stores 2.18 times the fine operator where the geometric one stores 1.48. What it is is applicable without a human deciding anything, on a problem class where the deciding cannot be done once for the whole grid.

What is asserted, and what refuses #

It rediscovers geometric coarsening: 31 of 63, none adjacent, weights ½ 1 ½ to exactly zero error, zero promotions.

Interpolation reproduces a constant on every interior row of two matrices, to 10⁻¹⁰.

The coarse operator stays symmetric to 10⁻¹⁴ and still annihilates a constant to 10⁻¹².

It discovers semi-coarsening: zero strong x-couplings, 31 of 31 uniform rows, 15 kept, none adjacent.

The threshold switches at ε, with every threshold above it giving a semi-coarsening and every one below giving a full one, and the fast side asserted to be fast.

And the sign matters, refused. A two-by-two matrix with one large positive off-diagonal is fed to strength(), and the assertion that it finds a strong connection is fed to rejects(). A strength measure that ignored the sign would pass it, and would then interpolate the rotated operator with weights of the wrong sign.

A coarsening that keeps everything is refused too. With θ above one, nothing is strong, no point is dropped, and what would run is a dense direct solve reporting a convergence factor of zero — the most flattering number available and the most misleading. The solve throws instead.

The cost of asking #

One thing the construction is not: free to build. Everything above happens before the first cycle, and it is a real amount of work on a matrix that has not yet been solved once.

The strength sets are one pass over the entries. The first coarsening pass is a greedy selection with a priority measure that has to be maintained as points are decided; the second is a pass over every strongly connected F-pair. The interpolation is another pass over the entries. And then the triple product PᵀAP, which is the expensive one, because forming it touches every path of length two through the fine matrix restricted to the coarse set — the same mechanism that makes the coarse operators dense in the first place.

All of that repeats at every level. A geometric hierarchy needs none of it: its transfer operators are three numbers known in advance and its coarse operators can be written down from the grid.

So the comparison between the two methods has three columns rather than two — the convergence factor, the cost of a cycle, and the cost of the setup — and only the first two are measured anywhere in this phase. That matters most for a matrix solved once, where the setup cannot be amortised, and least for one solved repeatedly with the same structure, where the hierarchy is built once and reused. Which of those an application is in decides the answer, and this phase does not measure either.

What is left #

The cost, which is the next essay and is substantial: the coarse operators fill in, and on an unstructured matrix one level of the hierarchy comes out entirely dense.

Smoothed aggregation, the other main family of algebraic methods, which builds coarse degrees of freedom by grouping fine ones rather than by selecting a subset, and is not touched here.

The rotated operator, where this construction manages 0.844 and the geometric line smoother manages 0.784. That is the phase’s honest stopping point and the next essay says why: the matrix does not contain the anisotropy the equation has.