Aggregating what the matrix calls strong

The depth phase ended the anisotropy field on a failure and a diagnosis. Four methods on the 45°-rotated anisotropic operator — a point smoother, a line smoother, semi-coarsening, and classical algebraic multigrid — returning convergence factors of 0.784, 0.883 and 0.844, none of them usable. And the reason, measured before any solver ran: the standard nine-point discretisation of the rotated operator has axis couplings of −0.5005 and a coupling along the strong direction of −0.2498, so the matrix does not contain the anisotropy the equation has.

That diagnosis makes a prediction, and the prediction is the reason this essay exists. If the information was lost in the discretisation, then no better hierarchy can recover it, because every algebraic method reads the matrix. Smoothed aggregation is the standard answer to anisotropy — the method a practitioner reaches for after classical AMG — so it is the sharpest available test of a claim about where the failure lives.

It fails, at 0.789, on the operator where it returns 0.193 unrotated.

The method, in four steps, and only one of them matters #

Strength. The same graph classical AMG builds: j is strong for i when −aᵢⱼ is at least θ times the largest off-diagonal magnitude in row i. Reused rather than rewritten, so the two methods disagree about the hierarchy and not about the matrix.

Aggregation. Cover the points with disjoint aggregates. Vaněk’s two passes: seed an aggregate at any point whose whole strong neighbourhood is still free, then hand each leftover point to the neighbouring aggregate it is most strongly attached to. Where classical AMG chooses a subset of the points to keep, this chooses a partition of them — which is the difference the two names are about. That the result is a partition is checked rather than assumed: every point appears in exactly one aggregate on all three operators, and the tentative prolongator’s columns are orthonormal to 10⁻¹².

A tentative prolongator. One column per aggregate, holding a normalised constant. A constant is what these operators annihilate, and an interpolation that cannot represent it cannot represent the error a smoother leaves behind.

And one weighted Jacobi sweep applied to the columns of that prolongator, with ω = 4/(3ρ(D⁻¹A)). This is the method.

The sweep is the method #

Same aggregates, same cycle, same smoother, same number of levels. One weighted Jacobi sweep applied to each column of P, and a factor of 2.3 in the convergence rate.

What the sweep does is widen the columns. A piecewise constant that is one on an aggregate and zero on its neighbour becomes a function that decays across the boundary — so the coarse space contains something with a gradient in it, and the error a smoother leaves behind has gradients. Unsmoothed aggregation asks a piecewise-constant space to represent a smooth error, and 0.813 is what that costs.

The constant is worth a sentence because this site has met it from the other side. ω = 4/(3ρ) is Vaněk’s choice and it damps the top two thirds of the spectrum; on the one-dimensional Laplacian, ρ(D⁻¹A) = 2 and 4/(3ρ) is exactly ⅔ — which is the smoother weight the multigrid field measured as optimal, where its optimum is exactly ⅓ and a 2,000-point search misses it and reports 0.3335.

What the aggregates know #

With ε = 10⁻³ the x-couplings are 10⁻³ and the y-couplings are 1, so the strength test sees the y direction and nothing else. Measured on a 31×31 grid: the aggregates’ mean extent is 2.82 points in y and 1.00 in x, over 341 aggregates.

That is semi-coarsening — every grid column kept whole or dropped whole — derived from the entries of the matrix, by a method with no coordinate in it anywhere. The depth phase found the same thing for classical AMG and recorded it as the phase’s organising sentence: the algebraic method rediscovers what the geometric one was told. It survives the change of method, which is what makes it a property of the strength graph rather than of one coarsening algorithm.

And the change from columns to blobs is a cliff at the strength threshold, not a gradient.

Along the whole slider the aggregation takes exactly two values. At ε = 1, 0.5, 0.3 and 0.25 it is 43 aggregates at 2.81 × 2.81; at ε = 0.2, 0.1 and 10⁻³ it is 75 at 3.00 × 1.00. Nothing happens in between and nothing happens outside, which is what a threshold test produces: the strength graph is a yes-or-no about each edge, so a continuous change in the operator moves the hierarchy only when it carries an edge across θ.

That has a practical reading the essay’s own diagnosis needs. A method whose coarsening is a step function of the operator cannot degrade gracefully as a problem becomes mildly anisotropic — it does the isotropic thing until ε crosses θ and the semi-coarsened thing afterwards, and at ε = 0.3 it is building an isotropic hierarchy for an operator whose couplings differ by a factor of three.

The coarsening ratio is the other thing two grid sizes settle. Anisotropic, it is 225/75 = 3.00 and 441/147 = 3.00 — exactly three at both, which is semi-coarsening by three in y and by nothing in x. Isotropic, it is 5.23 and 5.38. So the method’s coarsening rate is a property of the strength graph and not of the grid it happens to be run on, which is the sentence the whole algebraic programme rests on and is checkable here in one division.

Which direction counts as strong is read from coefficients(ε, 0) rather than assumed, and that is not pedantry: the two-dimensional operator’s coefficients are a = ε along x and c = 1 along y, so ε ≪ 1 makes y the strong direction, and the first version of this measurement asserted the extents the other way round and failed. A measurement of the right number attached to the wrong axis is the easiest mistake in this field to make and the hardest to see afterwards.

And the interpolation carries a constant, in the interior #

P₀ carries a constant exactly by construction — the coarse vector whose entries are √(aggregate size) maps to the all-ones vector, checked to 1.1·10⁻¹⁶.

The smoothed P does not, quite, and the measurement is worth reporting rather than glossing: (I − ωD⁻¹A) applied to a constant leaves it alone wherever the row sum of A is zero, which is every interior row, and changes it where the row sum is not — the boundary. Measured: the deviation from one is 1.1·10⁻¹⁶ over the interior points and 0.34 at the worst boundary point.

So the near-null-space property the method rests on is exact where the operator is homogeneous and broken at the edge, by an amount that is a boundary effect. That is the same shape as the finding two essays away in the structure field, where a clustering result thins because its exception set is the boundary of the grid — and it is worth noticing that two unrelated methods have their one-dimensional exactness spoiled by the same geometric fact.

The rotation, and the prediction it tests #

The rotated operator’s coefficients are a = c = (ε + 1)/2 = 0.5005 and b = (ε − 1)/2 = −0.4995, so the stencil has four axis entries of −0.5005 and four diagonal entries of ∓0.2498. Two consequences, and both of them break the strength test rather than the solver:

The largest off-diagonal in every interior row is an axis coupling. The strength test compares −aᵢⱼ against the largest in the row, so it marks the axis directions strong and the diagonal directions weak — and the diagonal directions are where the anisotropy runs.

And two of the four diagonal entries are positive. Every classical strength measure ignores positive off-diagonals entirely, for a good reason amg.js states: for a matrix with zero row sums, a positive off-diagonal is a coupling that pushes two values apart rather than together, and treating it as strong produces interpolation weights of the wrong sign. So half of the couplings along the strong direction are not merely under-weighted, they are invisible.

A method reading that matrix has no way to find the direction, because the direction is not in the matrix. Smoothed aggregation reads that matrix. It returns 0.789.

Where the failure actually begins #

Two angles is enough to establish that the method fails and not enough to say when. Sweeping the angle, k = 31, ε = 10⁻³:

From 0° to 25° the hierarchy does not change at all. The same 341 aggregates, the same 2.82 points tall and 1.00 wide, the same semi-coarsening — and the convergence factor triples over that range, with most of the loss arriving by ten degrees.

Over the whole of that range the strength test is right. It reports y as the strong direction, y is the nearest grid direction to the anisotropy, the diagonal couplings are 0.088 of the axis ones at 10° — well under θ = 0.25 — and the aggregates follow the direction the test found. Nothing is being fooled. The method degrades because the anisotropy points ten degrees away from the only direction the grid has, and there is no aggregate shape that can point ten degrees. A line smoother has the same limitation for the same reason, and smoothing a whole line at once is the essay that prices it.

The 45° diagnosis is about the second regime, and θ is the boundary. The diagonal coupling relative to the largest axis coupling crosses 0.25 between 25° and 30°, and the aggregation changes abruptly there rather than gradually: 341 columns become 121 blobs, one point wide becomes 2.82 points wide. Semi-coarsening is abandoned in a single step. That is the same step change the aggregation shows as a function of ε, appearing here as a function of the angle, and it is what the stencil argument above describes.

Separating the two regimes matters because they point at different repairs, and the essay’s recommendation only reaches one of them. A stencil aligned with the anisotropy addresses the second: at 45° the discretisation genuinely loses the direction, and putting it back is a discretisation question. It cannot address the first. At 10° the operator is nearly separable, the couplings are faithful, the direction is in the matrix and correctly found — and the method has already lost sixty per cent of its rate. Only a mesh that follows the flow reaches that regime, and that regime is where most of the degradation happens. That is the same conclusion exact along one axis reaches from the convection field: a scheme whose competence is aligned with the grid loses it the moment the problem is not.

One irregularity worth recording rather than smoothing: the factor is not monotone in the angle. It rises to 0.636 at 15°, sits at 0.634 at 20°, and falls to 0.542 at 25° before the step. Whatever governs that is finer than the direction argument, and the direction argument is what this essay claims.

What a failure that was predicted is worth #

Two things, and it is worth separating them from the disappointment.

The diagnosis is confirmed rather than repeated. The depth phase’s claim was that the failure at 45° belongs to the discretisation. Its evidence was that four methods failed and that the stencil’s couplings pointed the wrong way. Adding a fifth failing method is weak evidence on its own; adding the one method whose entire purpose is anisotropy, which succeeds on the same anisotropy when it is axis-aligned, is much stronger. The variable that changed is the rotation and not the method.

And the fix is now locatable. If the information is missing from the matrix, the repair has to put it back before any solver sees it — which means the discretisation. A nine-point stencil aligned with the anisotropy rather than with the grid, or a mesh that follows the flow, or a stabilisation that knows the direction, all do that; a better hierarchy cannot. That is a conclusion about where to spend effort, and it is the kind of conclusion a failure is good for.

This is the same move the previous essay makes in the convection field, and the two are worth reading together: there, a scheme succeeds spectacularly because the answer’s shape was built into it, and fails on the problem next door. Here, a method fails because the operator’s shape was lost before it arrived. Both are the depth phase’s sentence — where a method’s competence was put into it by hand — measured from the two sides.

What it costs when it works #

The honest accounting, since a method that converges at 0.193 is not free:

Operator complexity is the total nonzeros across the hierarchy divided by the fine matrix’s, and it is the quantity that decides whether an optimal-complexity method is actually cheap. The smoothing widens P, which densifies every coarse operator below it, and 1.34 is a good number — classical AMG on the same grid runs higher, and the depth phase measured 17.7 on a graph Laplacian, with one level of forty-one unknowns entirely dense.

So this method is cheaper than classical AMG per cycle and converges more slowly per cycle, which is the standard trade between them and is measured here rather than repeated.

What the two passes each do, in counts #

The aggregation is two passes and the split between them is a measurement rather than an implementation detail — a method whose second pass is doing most of the work is a method whose first pass found little structure.

Read the last row against the first two. On the rotated operator the first pass opens thirty aggregates and the second pass has to place forty-five leftover points — half again as many as were seeded — where on the aligned operator it places fifteen. The strength graph the first pass walks is finding fewer clean neighbourhoods, because the couplings it can see are the axis ones and those do not agree about which points belong together.

That is a diagnostic available before any solve, and it is the second of the two the stencil supplies. A method whose aggregation is mostly second-pass adoption is a method operating on a graph that does not describe the problem.

The near-null space, and what “carries a constant” means precisely #

The tentative prolongator carries a constant exactly: the coarse vector whose entries are √(aggregate size) maps to the all-ones vector, to 1.1·10⁻¹⁶. That is by construction — the columns are normalised indicators of a partition — and it is why the partition property is checked rather than assumed.

The smoothed prolongator does not, quite, and the measurement is worth stating rather than glossing. (I − ωD⁻¹A) applied to a constant leaves it alone wherever the row sum of A is zero, which is every interior row of these operators, and changes it where the row sum is not — the boundary. Measured on the aligned operator: the deviation from one is 1.1·10⁻¹⁶ over the interior and 0.34 at the worst boundary point.

So the property the method’s theory rests on is exact in the interior and broken at the edge by a bounded amount. Both halves are asserted, and the reason to assert the second is that it is the difference between a claim that is true of an infinite grid and a claim that is true of the matrix in front of the method.

That is the same shape as the structure field’s loss two fields away, where a clustering result weakens in two dimensions because its exception set is the boundary of the grid rather than the ends of a vector. Two unrelated methods, the same geometric fact.

What the aggregation would have to see #

If the diagnosis is right — and the essay above is an experiment designed to test it rather than an argument for it — then the repair belongs in one of two places, and neither is the hierarchy.

A strength measure that sees positive couplings. Every classical measure compares −aᵢⱼ against the largest in the row and ignores positive off-diagonals, for a reason amg.js states: for a matrix with zero row sums a positive off-diagonal pushes two values apart. The rotated stencil has two positive diagonal entries per interior row, so half of the couplings along the anisotropy are not under-weighted but invisible. Affinity and algebraic-distance measures do see them, and whether they recover the 45° direction on this stencil is a measurement rather than an argument.

Or a discretisation that keeps the direction. The nine-point stencil used here is the standard one and it is where the information is lost. A stencil aligned with the anisotropy, or a mesh that follows it, would hand every method in this field a matrix that contains what the equation contains — and the measurement that matters is whether the four failing methods then return to the 0.19 they reach when the anisotropy happens to lie along an axis.

What is left #

Energy-minimising interpolation, the other half of the depth phase’s deferral, which chooses P by minimising the energy of its columns subject to reproducing the near-null space rather than by applying one Jacobi sweep. It is the natural next step from the smoothing sweep above — a sweep is one iteration of exactly that minimisation — and it is the method most likely to do better on the rotated operator, because it is the least dependent on the strength graph.

A strength measure that sees positive couplings. The affinity and algebraic-distance measures do, and whether they find the 45° direction in this stencil is a measurement rather than an argument. If they do, this essay’s conclusion needs narrowing to classical strength measures; if they do not, the diagnosis is stronger than stated.

And the discretisation itself, which is where the fix has to go and which this whole field has so far treated as given. The nine-point stencil used everywhere here is the standard one, and the measurement that matters is whether a stencil built to align with the anisotropy restores any of the four methods — including this one — to the 0.19 they reach when the anisotropy happens to lie along an axis.

What calling the matrix requires #

Deciding which connections are strong means comparing entries within a row. That is the operation an operator supplied as a subroutine cannot perform, which is why every method in this field needs an assembled matrix.