Every essay — page 4
Iterating, instead of factorising
At scale nobody factorises, and the reason is not accuracy — it is that the factors of a sparse matrix are not sparse. What replaces elimination is a sequence of approximations, and the question changes shape: not what the residual of the factorisation is, but how fast the error falls and whether the thing you can measure tracks the thing you cannot. On the model problem every rate is known in closed form before anything runs.
The error smoothing cannot reach
One weighted Jacobi sweep multiplies every mode of the error by a number, and the number is a sine. Half the modes are cut by three or better, and the other half come back at 0.999 — which is not a failure of the method but the fact the whole of multigrid is built on.
The same problem on a coarser grid
Restriction, the coarse operator and interpolation are three matrices with nine distinct entries between them. Two of the three are each other's transpose, and their product with the fine operator is the coarse discretisation exactly — not approximately, entry for entry, at every level.
A rate that does not notice the size
The V-cycle reduces the residual by a factor of ten a cycle at fifteen points and at a hundred and twenty-seven. Jacobi on the same four problems goes from 0.981 to 0.9978, climbing towards one. One of those is a constant and the other is an exponent, and that is the whole distinction the field turns on.
The coarse problem is a different problem
In one dimension the Galerkin coarse operator is the coarse discretisation, entry for entry — this site asserted it. In two dimensions a five-point operator produces a nine-point coarse one, so the recursion solves a different discretisation at every level below the first, and converges at 0.20 a cycle regardless.
A direction the smoother cannot see
Give the Laplacian a strong direction and multigrid stops working — from 0.2016 a cycle to 0.9565 — with every component unchanged and the condition number identical to twelve digits. The problem did not get harder. The link between the method's two halves broke.
Smoothing a whole line at once
Solve every grid line in the strong direction exactly rather than sweeping over it, and the smoothing factor goes from 0.9993 back to 0.3340 — which is the one-dimensional answer, on a problem that is not one-dimensional. The repair replaces one ε in the closed form by a one.
Coarsening in one direction only
Leave the smoother alone and halve only the strong direction, and the smoothing factor is 0.3340 — identical to line relaxation's, at every anisotropy and every weight, to twelve digits. The convergence factors are then a factor of three apart, and at 45° both repairs fail outright.
The coarse grid the matrix chooses
Given a tridiagonal matrix and no information about a grid, the coarsening keeps every other point and derives the weights ½, 1, ½ — the operators the geometric method was handed. Given the anisotropic operator, it discovers semi-coarsening, in the right direction, without a coordinate.
A hierarchy with no grid behind it
On a graph Laplacian the algebraic V-cycle converges at 0.199 a cycle, its grid complexity is an unremarkable 3.05, and its operator complexity is 17.7 — one level of forty-one unknowns is entirely dense. The number people quote is the one that does not measure the work.
The formula that was already optimal
Ask for the interpolation that minimises the energy of its own columns and the answer is the classical AMG formula — to zero at every row of the one-dimensional Laplacian, and to four digits in two dimensions. On the operator rotated to 45° the two part company, and the gap between them is a diagnostic that needs no reference solution.
The stencil that is not symmetric
Past a cell Péclet number of exactly one — measured by bisection at 1.0000000000000002 — the central-difference solution of a convection–diffusion problem oscillates from point to point and leaves the interval the equation guarantees, at 16 of 31 grid points. It is the exact solution of its own linear system, to 4.6·10⁻¹⁸. No solver was involved.
The diffusion that makes the answer exact
Upwinding adds h/2 of artificial diffusion. Central differencing adds none. Add ε·ξ·Pe with ξ = coth(Pe) − 1/Pe and the computed solution is the exact one at every grid point, to 2.4·10⁻¹⁷ — at every Péclet number, on the problem it was derived from and on no other.
Aggregating what the matrix calls strong
The depth phase measured every method it had on the 45°-rotated anisotropic operator — 0.784, 0.883, 0.844 — and diagnosed the failure as being in the discretisation rather than in the hierarchy. Smoothed aggregation is the standard answer to anisotropy. It returns 0.789.
One sequence and two recurrences
CGLS and LSQR compute the same iterates — the minimiser over a space is unique, so there is nothing to choose between them in the algebra. At κ = 10⁶ they cost 42 steps and 47. At κ = 10¹⁰ they cost 110 and 209, across four seeds, and the quantity that separates them is the orthogonality of a basis neither of them keeps.
Exact along one axis
The tuned diffusion makes the answer exact at every node, and in two dimensions it holds at exactly one flow angle. Five degrees off the grid the relative error goes from 1.2·10⁻¹⁴ to 6.9, and by twenty degrees the scheme is worse than the upwinding it was built to improve on.
The direction the diffusion does not go
Streamline diffusion adds τbbᵀ, a rank-one tensor that annihilates every direction across the flow. That is the design. The price is 18, 50 and 48 nodes where the computed solution leaves the interval the equation guarantees — and half a coefficient of crosswind diffusion halves the error at every angle while costing exactly nothing where the scheme was exact.
An operator with no entries
At the sizes where linear algebra is expensive the matrix does not exist. What exists is a subroutine that returns Av. Every Krylov method survives that unchanged; every algorithm that reads an entry disappears. And the derivative such a code computes is accurate to ten digits instead of sixteen, which turns out to cost nothing at all.
Sparsity, and what elimination costs
Eliminating a variable couples everything it touched to everything else it touched, and every one of those couplings is an entry that was zero in the matrix and is not zero in its factor. On the same matrix one elimination order gives a factor of a thousand entries and another gives ten thousand — the two factorisations equally accurate, one of them fitting in memory. Nothing numerical chooses between them.
The factor is not sparse
A sparse matrix has a factor that is not sparse, and the gap between them is the entire reason iterative methods exist. The entries elimination creates can be counted before any arithmetic runs, from the graph alone.
The order decides the memory
Four elimination orderings on one matrix give factors of 1,739, 1,354, 1,413 and 1,026 entries. All four factorisations are exact, all four return the same answer, and the one with the better asymptotics is not the one that wins.
Two ends of the same arrow
One matrix, one row moved from the front of the elimination order to the back, and the factor goes from completely dense to no fill at all. Both factorisations are exact to rounding, and nothing numerical chose between them.
Structure and stability stop being separable
The sparsest variable to eliminate on this matrix has a diagonal entry of 10⁻¹². Eliminating it produces the smaller factor, reproduces the matrix to 3.8·10⁻¹⁷ — better than pivoting does — and returns an answer wrong in the fifth digit.
A threshold between fill and growth
One number decides how small a pivot an elimination will accept. At 0.001 the factor holds 172 entries and the matrix grows by 1,330; at 1 it holds 260 and grows by 1.2. The libraries ship 0.1, and the measurement says why.
What the symbolic phase can only bound
Without pivoting, the fill can be computed from the graph and the count is exact — 233 predicted, 233 measured. With pivoting it is 233 predicted and 242 measured, and what survives is a bound that is right at every threshold and loose by 1.7 times at the largest grid drawn.
Structure, and the solver that cannot see it
The sparsity field is about a matrix most of whose entries are zero. This is the other way a matrix can be small: a circulant is n numbers and a Toeplitz matrix is 2n − 1, and neither of them has a single zero entry, so nothing in the sparsity field applies. What the structure buys is an exact spectrum, a condition number in closed form and a solve in n log n — and what it does not buy is a well-conditioned problem, which is the half that gets summarised away.