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ᵢⱼ  ≥  θ · max over k ≠ i of (−aᵢₖ)

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 two frames so far are the two ends. What the slider settles is whether there is anything in between them.

Two hundred times the anisotropy and not one number has moved. The next step is a tenth of that change and moves all four:

Ten stops and two answers. Every reading at ε ≤ 0.2 is identical and every reading at ε ≥ 0.3 is identical, to the last integer — the same 0, 110, 55, 11 six times and the same 110, 110, 61, 0 four times. There is no intermediate regime at all, on a slider spanning three orders of magnitude.

The step is at θ, exactly. The threshold is 0.25, and the coupling ratio is ε; ε = 0.2 is below it and ε = 0.3 is above it, and that is the whole of what happens. The coarse grid the matrix chooses is not a function of the anisotropy — it is a function of one comparison between the anisotropy and a number somebody typed into a parameter file.

And all 110 couplings switch together. Not a few, not a growing fraction: 0 to 110 in one step, because on a uniform grid every x-coupling has the same ratio and a single threshold either admits all of them or none. So the parameter has no interior at all here, and a code tuning θ on this problem is not adjusting a balance — it is choosing between semi-coarsening and full coarsening with nothing between.

The two grids differ by six points, 55 against 61 of 121, which is worth noting beside the rest: the strength graph changes completely and the coarse grid it produces changes by five per cent. What changes is the shape — 11 of 11 rows treated whole below the step and 0 of 11 above it — and the count barely moves.

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, which is the arrangement the same problem on a coarser grid describes and the arrangement this construction is doing without.

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. The first two are measured throughout this collection and the third was, for a long time, left as a remark. It is worth having the number, because the remark understates it.

Timing the two halves against each other on the same operator, tolerance 10⁻¹⁰, each repeated after a warm-up:

For a matrix solved once, more than half the work happens before the first residual is computed. That is not the impression “there is also a setup cost” leaves, and it is the impression the sentence was doing duty for.

The second column of that table is the one that decides whether it matters, and it is flat. Across a factor of 6.8 in the number of unknowns the ratio stays near ten on the isotropic operator and near seven on the anisotropic one. So the setup is a constant multiplier on the solve, not an overhead that a larger problem dilutes — the two halves have the same complexity, which is what one would hope for and is not the same as one of them being negligible. A code that plans to amortise the setup has to amortise it over repeated solves, not over a bigger grid.

Wall-clock is a poor instrument for an absolute cost and a fair one for this ratio, since both halves run in one process against one representation of the matrix; the assertion’s bands are wide for that reason. But there is a deterministic quantity underneath it that says the same thing more sharply, and it holds a surprise. Operator complexity — the nonzeros of the whole hierarchy over the nonzeros of A, which is exactly what the setup builds and stores — behaves differently on the two operators:

The isotropic column is a constant to three digits over a factor of seven in n, which is the textbook statement. The anisotropic column climbs by 20% over the same range, and the depth goes 5, 6, 7, 7 against 4, 5, 5, 6. Semi-coarsening is the reason: it halves in one direction only, so each level removes half as much as a full coarsening does and more levels are needed to reach the same coarse size, and every one of them adds its operator to the total. That is coarsening in one direction only seen from the storage side rather than the convergence side.

That is the cost of the right decision, and it is worth naming as such. The algebraic construction finds semi-coarsening without being told to, which is the result this essay is about; semi-coarsening is what makes the convergence factor an order of magnitude better; and semi-coarsening is also what makes the hierarchy grow. The method is not paying for a mistake here. It is paying for the anisotropy, in the one place the convergence factor does not show it — which is the same shape as a block size is a property of the machine, where the quantity that decides the answer is not the one the analysis is written in.

assertTheSetupIsAboutTenCycles in the algebraic library measures both halves and both tables, and requires the flat column to stay flat and the climbing one to climb.

Which regime an application is in still decides the answer: the setup matters most for a matrix solved once and least for one solved repeatedly with the same structure. What has changed is that “most” is now a number, and it is more than half.

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.

What choosing from the matrix costs #

Choosing a coarse grid from the matrix entries is what makes this method algebraic, and it is what makes it unavailable where there are no entries. The trade is worth measuring on something simpler.