Algebraic multigrid preconditioning

A worked comparison of three ways to solve the same discrete Poisson system: a direct sparse factorization, plain (unpreconditioned) conjugate gradient, and conjugate gradient preconditioned by algebraic multigrid (AMG). Every number below was produced by the code shown.

Problem

\[-\Delta u = g \text{ in } \Omega, \qquad u = u_{\text{exact}} \text{ on } \partial\Omega, \qquad \Omega = (0,1)^2\]

with the manufactured solution $u_{\text{exact}}(x, y) = e^{x+y}$, so $g = -2 u_{\text{exact}}$. This is deliberately not a trigonometric solution such as $\sin(\pi x)\sin(\pi y)$: on a uniform grid that happens to be (very nearly) a single eigenmode of the discrete Laplacian, which both plain and preconditioned CG then solve in a handful of iterations regardless of mesh size – a measurement made once already, that looked like a working comparison and said nothing about the preconditioner at all. $e^{x+y}$ excites the discrete spectrum broadly, so the iteration counts below reflect the operator's actual conditioning.

Why this needs a preconditioner at all

The condition number of the assembled Laplacian scales as $\mathcal{O}(h^{-2})$, so an unpreconditioned Krylov method needs $\mathcal{O}(h^{-1})$ iterations – doubling, very roughly, every time the mesh is refined by a factor of two. Algebraic multigrid builds a coarse-grid hierarchy directly from the graph of the assembled matrix and gives a Krylov method grid-independent, $\mathcal{O}(1)$ iterations instead. amg_preconditioner wraps AlgebraicMultigrid.jl's two constructions, :smoothed_aggregation (the default) and :ruge_stuben.

Assembling a symmetric system

AMG is built for symmetric positive-definite matrices. assemble(a::BilinearForm; dirichlet = ...) alone does not produce one: a Dirichlet row becomes the identity, but the matching column is left alone, so the assembled matrix is not exactly symmetric even though the underlying operator is. Passing symmetrize = true to the two-form assemble (or to linear_problem/solve below) restores that symmetry by folding the removed columns into the right-hand side – do this before handing the matrix to AMG, not after.

using Bramble
using LinearSolve
using AlgebraicMultigrid

uex(x) = exp(x[1] + x[2])
rhs(x) = -2 * uex(x)

function poisson_system(n)
    Ω = domain(interval(0.0, 1.0) × interval(0.0, 1.0))
    Ωₕ = mesh(Ω, (n, n), (true, true))
    Wₕ = gridspace(Ωₕ)
    bcs = dirichlet_constraints(Ω, :boundary => uex)

    a = form(Wₕ, Wₕ, (u, v) -> inner₊(∇ₕ(u), ∇ₕ(v)))
    gₕ = element(Wₕ)
    avgₕ!(gₕ, rhs)
    l = form(Wₕ, v -> innerₕ(gₕ, v))

    A, F = assemble(a, l; dirichlet = bcs, symmetrize = true)
    return a, l, bcs, A, F
end

a, l, bcs, A, F = poisson_system(65)
issymmetric(A)

true

Three ways to solve it

Direct: A \ F, a sparse LU/Cholesky factorization – exact up to round-off, and perfectly fine at this size.

u_direct = A \ F

4225-element Vector{Float64}:
 1.0
 1.0157477085866857
 1.0317434074991028
 1.0479910020166328
 1.0644944589178593
 1.0812578074490395
 1.0982851403078258
 1.1155806146424807
 1.1331484530668263
 1.1509929446911764
 ⋮
 6.5208191203301125
 6.623507079583559
 6.727812138894691
 6.833759763883972
 6.941375821197036
 7.050686584819912
 7.161718742493711
 7.274499402230307
 7.38905609893065

Plain CG: no preconditioner.

sol_cg = solve(LinearProblem(A, F), KrylovJL_CG(); reltol = 1e-8, abstol = 1e-10)
sol_cg.iters

AMG-preconditioned CG: build the hierarchy once with amg_preconditioner, wrap it with AlgebraicMultigrid.aspreconditioner, and pass it as Pl.

P = aspreconditioner(amg_preconditioner(A))
sol_amg = solve(LinearProblem(A, F), KrylovJL_CG(); Pl = P, reltol = 1e-8, abstol = 1e-10)
sol_amg.iters

The same three lines collapse into one call through linear_problem's companion solve(a::BilinearForm, l::LinearForm; ...), which assembles, solves and unwraps the result to a VectorElement directly – preconditioner = :amg reaches amg_preconditioner internally, so there is nothing to build by hand:

uₕ = solve(a, l; dirichlet = bcs, symmetrize = true, solver = KrylovJL_CG(), preconditioner = :amg)
normₕ(uₕ .- Rₕ(space(uₕ), uex))

4.308760619495785e-8

The point of the comparison: plain CG's iteration count against AMG's, on the same manufactured problem, as the mesh is refined.

function iters(n)
    _, _, _, A_n, F_n = poisson_system(n)
    prob = LinearProblem(A_n, F_n)

    plain = solve(prob, KrylovJL_CG(); reltol = 1e-8, abstol = 1e-10).iters

    P_n = aspreconditioner(amg_preconditioner(A_n))
    amg = solve(prob, KrylovJL_CG(); Pl = P_n, reltol = 1e-8, abstol = 1e-10).iters

    return plain, amg
end

ns = (16, 32, 64, 128)
results = iters.(ns)
plain_iters = first.(results)
amg_iters = last.(results)

for (n, p, m) in zip(ns, plain_iters, amg_iters)
    println("n = $(lpad(n, 3))   plain CG = $(lpad(p, 4))   AMG-CG = $(lpad(m, 3))")
end

n =  16   plain CG =   45   AMG-CG =   6
n =  32   plain CG =   94   AMG-CG =   7
n =  64   plain CG =  188   AMG-CG =   7
n = 128   plain CG =  371   AMG-CG =   9

plain CG roughly doubles at each refinement – the O(h^-1) growth the condition number predicts – while AMG-CG stays within a handful of iterations across a 64-fold increase in degrees of freedom.

Algebraic multigrid preconditioning

Problem

Why this needs a preconditioner at all

Assembling a symmetric system

Three ways to solve it

Iteration counts under refinement