Multigrid Methods

• Multigrid methods are iterative techniques for solving large linear systems derived from elliptic PDEs, focusing on reducing high and low-frequency errors through fine and coarse-grid corrections, ensuring efficient convergence.
• Multigrid gives accurate solutions for varied problems like finite-element formulations, adaptive meshes, sodium-based simulations, lattice fermions, and uncertainty-sampled saddle-point systems by optimizing smoothing of error representation at multiple scales.
• Applications, such as adaptive meshes, spectral transformations, accelerators, and multi-grid problem-solving helps tackle problems like uncertainties, approximate Dirac operator in Grid computing, and thickened assemblies.

Multigrid is a family of iterative and preconditioning methods for large linear systems, especially those obtained by discretizing elliptic partial differential equations. Its defining mechanism is the complementary use of relaxation and a hierarchy of coarser representations: relaxation rapidly reduces oscillatory, high-frequency error, whereas coarse-grid correction represents and removes smooth, low-frequency error. In a variational setting for an SPD system $Au=f$, the solution minimizes $E(u)=\frac12u^TAu-f^Tu$, and multigrid combines inexpensive approximate energy reduction on the fine grid with recursively computed corrections on coarser grids. Geometric multigrid constructs levels from spatial meshes; algebraic multigrid derives them from matrix structure; related methods extend the concept to unstructured surfaces, adaptive meshes, lattice fermions, uncertainty-sampled saddle-point systems, modal representations, and learned transfer or smoothing operators.

1. Fundamental formulation and error decomposition

For an elliptic problem

$Lu=f,$

discretization on a grid of mesh size $h$ produces

$L^h u^h=f^h.$

When $L^h$ is symmetric positive definite,

$(L^h)^T=L^h,\qquad x^TL^hx>0$

for every nonzero $x$. The system is equivalently the minimization of a strictly convex quadratic energy. If $u_\ast=A^{-1}f$ and $e=u-u_\ast$, then

$E(u)=\frac12u^TAu-f^Tu$0

so the energy norm is

$E(u)=\frac12u^TAu-f^Tu$1

The residual is the energy gradient. Depending on convention, it is written either as $E(u)=\frac12u^TAu-f^Tu$2 or as $E(u)=\frac12u^TAu-f^Tu$3. With the latter convention, the correction equation is

$E(u)=\frac12u^TAu-f^Tu$4

This distinction changes only signs if used consistently.

The central multigrid observation is spectral. For elliptic discretizations, oscillatory error generally corresponds to large eigenvalues, while slowly varying error corresponds to small eigenvalues. A relaxation method therefore reduces oscillatory error efficiently but acts slowly on smooth error. For example, Richardson relaxation has the form

$E(u)=\frac12u^TAu-f^Tu$5

and its error-propagation matrix is

$E(u)=\frac12u^TAu-f^Tu$6

With $E(u)=\frac12u^TAu-f^Tu$7, an eigencomponent associated with eigenvalue $E(u)=\frac12u^TAu-f^Tu$8 is multiplied by

$E(u)=\frac12u^TAu-f^Tu$9

Large-eigenvalue components are rapidly damped; small-eigenvalue components remain. Algebraic smoothness is therefore more general than geometric smoothness: an error is algebraically smooth when applying the operator produces a relatively small residual compared with the error itself.

2. Grids, transfer operators, and multigrid cycles

A two-grid method supplements fine-grid relaxation with coarse-grid correction. Let $Lu=f,$0 prolongate a coarse vector to the fine grid and $Lu=f,$1 restrict a fine residual to the coarse grid. A typical cycle consists of pre-smoothing, residual computation, restriction, approximate coarse solution, prolongation, correction, and post-smoothing:

$Lu=f,$2

$Lu=f,$3

$Lu=f,$4

$Lu=f,$5

For variational methods, restriction is often chosen as the transpose of prolongation,

$Lu=f,$6

and the coarse operator is formed by the Galerkin product

$Lu=f,$7

The Galerkin correction minimizes the fine-grid energy over the coarse representable subspace. After correction, the remaining error satisfies the orthogonality condition

$Lu=f,$8

The principal recursive cycles are:

1. V-cycle: one recursive coarse-grid visit at each level.
2. W-cycle: two recursive visits, providing stronger coarse correction at greater cost.
3. F-cycle: an initialization-oriented cycle between full multigrid and V-cycle behavior.
4. Full multigrid (FMG): begins on the coarsest grid, solves there, prolongates toward the finest grid, and applies cycles at successive levels.

FMG targets discretization-level accuracy rather than merely algebraic residual reduction. Its objective is to compute a solution whose error is comparable to the discretization error at a cost comparable to a few fine-grid matrix operations.

In finite-element formulations, transfer operators may follow exactly from basis-function relations. For ghost finite elements, this eliminates residual splitting between interior, Dirichlet-boundary, and Neumann-boundary equations. For unfitted domains, Galerkin transfer remains consistent with Nitsche boundary and penalty terms, provided stabilization is selected to preserve coercivity.

3. Geometric, algebraic, and semi-structured multigrid

Geometric multigrid uses explicitly constructed meshes such as

$Lu=f,$9

and rediscretizes or transfers the PDE between levels. It is highly efficient when geometry, discretization, and coarsening are regular. A geometric multigrid library for quadtree and octree adaptive meshes supports Cartesian and cylindrical geometries, second-order five- and seven-point operators, periodic, Dirichlet, Neumann, and three-dimensional free-space boundary conditions. Its AMR hierarchy is 2:1 balanced, and its matrix-free implementation has demonstrated strong scaling to 1792 cores (Teunissen et al., 2019).

Geometric coarsening becomes difficult on unstructured meshes, curved surfaces, locally refined finite-element grids, and domains represented implicitly by level sets. Auxiliary-space multigrid addresses this by constructing a geometrically organized auxiliary hierarchy while retaining a smoother on the original finite-element space. For shape-regular, locally refined unstructured meshes, cluster-tree construction, closure, storage, and multilevel assembly have complexity $h$0. The resulting method has convergence rate

$h$1

and condition-number estimate

$h$2

This is nearly optimal rather than uniformly optimal, but it applies without requiring the original mesh to be globally quasi-uniform (Grasedyck et al., 2014).

Algebraic multigrid (AMG) constructs coarse spaces from matrix and graph information. Aggregation AMG forms aggregates using a maximal independent set on a graph of strong connections. With one near-null-space vector, unsmoothed piecewise-constant interpolation has one nonzero per fine-grid row,

$h$3

This makes setup, storage, restriction, prolongation, and coarse-matrix construction inexpensive. Because unsmoothed aggregation alone may not provide grid-independent convergence, recursively applied coarse Krylov iterations, called K-cycles, can be used on upper levels, while V-cycles handle lower levels. A hybrid K-cycle/V-cycle preconditioner has been used with flexible GMRES and achieves nearly constant iteration counts for tested systems growing from 1.1 million to 13.5 million unknowns (Gandham et al., 2014).

Semi-structured AMG, or SSAMG, retains logically rectangular parts and represents only inter-part couplings in an unstructured component. If

$h$4

where $h$5 is the structured stencil component and $h$6 contains inter-part connections, then coarse operators retain the decomposition

$h$7

Semi-coarsening directions are selected independently for each part using matrix-derived metrics. This permits local adaptation to anisotropy while preserving regular stencil computation in the bulk. In reported tests, SSAMG setup was up to approximately four times faster than BoomerAMG and remained applicable to multi-part, anisotropic, and structured adaptive-mesh problems for which a single structured solver was inadequate (Magri et al., 2022).

4. Irregular geometry, adaptive discretization, and difficult smoothers

Multigrid convergence depends not only on the abstract hierarchy but also on the quality of transfer operators, smoothers, and intermediate representations.

For triangulated curved surfaces, ordinary barycentric interpolation fails because fine vertices generally do not lie in the plane of corresponding coarse triangles. Intrinsic surface multigrid constructs a bijective fine-to-coarse map through successive local parameterizations during edge collapses. Each fine vertex is interpolated from three coarse vertices, so the prolongation has three nonzeros per row. The geometry-only hierarchy can be reused for cotangent Laplacians, squared Laplacians, squared Hessians, smoothing systems, deformation, and simulation. Reported applications include speedups of $h$8 and $h$9 for squared-Hessian problems with 200,000 and 3 million vertices, respectively (Liu et al., 2021).

Gravo MG uses graph Voronoi coarsening rather than constructing watertight coarse meshes. Coarse points are selected by spatially uniform subsampling, graph Voronoi cells define coarse adjacency, and local candidate triangles provide barycentric interpolation. The method applies to triangle meshes, polygon meshes, point clouds, and nonmanifold meshes. Its prolongation has at most three nonzeros per row, and hierarchy construction is reported to be on average 36 times faster than the intrinsic multigrid construction of Liu et al. (Wiersma et al., 2023).

Adaptive and unfitted domains create different difficulties. In ghost finite elements, cut cells can produce poorly conditioned contributions and boundary-localized error. The stabilization parameter must be sufficiently large for coercivity but not unnecessarily large, because excessive Nitsche penalization damages smoothing. In two dimensions, additional smoothing restricted to cut cells is inexpensive because the number of cut cells grows like boundary measure divided by $L^h u^h=f^h.$0, whereas the total number of cells grows like $L^h u^h=f^h.$1. With local stabilization and four additional cut-cell smoothing iterations, reported convergence factors are approximately

$L^h u^h=f^h.$2

for several curved geometries (Dilip et al., 8 May 2025).

Low-quality cells create a related but distinct local failure. Standard smoothers may damp error throughout a mesh while leaving a concentrated residual near degenerate elements. A global–local smoother supplements global relaxation with exact local residual corrections on patches surrounding cells whose normalized radius ratio satisfies

$L^h u^h=f^h.$3

The combined operator has the symmetric form

$L^h u^h=f^h.$4

This restores high-quality-mesh convergence in two- and three-dimensional Poisson and elasticity experiments, with particularly strong effects for quadratic $L^h u^h=f^h.$5 elements (Chen et al., 2024).

For strongly anisotropic operators, pointwise smoothers may fail because error is strongly coupled along lines or planes. Alternating line Gauss–Seidel is therefore used for multipoint flux approximations of Darcy problems on rough quadrilateral grids. Local Fourier analysis shows that point Gauss–Seidel can fail for highly anisotropic tensors, whereas alternating line relaxation remains convergent. The method combines piecewise-constant prolongation, Wesseling/Khalil restriction, Galerkin coarsening, and line relaxation (Arrarás et al., 2024).

5. Multigrid in many-core and lattice-QCD computations

Multigrid algorithms have favorable algorithmic complexity but can map poorly to accelerators when coarse grids contain fewer points than processors. High-occupancy multigrid (HOMG) addresses this on GPUs by representing logical coarse operations over the entire fine-grid domain. Coarse-grid points perform relaxation, non-coarse points perform interpolation, and the two operations are fused. Data remain on the GPU, regular memory access is retained, and occupancy remains high.

HOMG uses damped Jacobi and a modified full multigrid cycle. At fixed accuracy, the reported work scales linearly with the number of grid nodes, and large-grid experiments approach the ideal scalability curve. The method was demonstrated on a structured two-dimensional Poisson problem, with reported errors of approximately 5% for $L^h u^h=f^h.$6 and 1.6% for $L^h u^h=f^h.$7. Its evidence is limited to structured two-dimensional single-precision experiments and does not establish universal superiority over conventional or hybrid implementations (Becerra-Sagredo et al., 2011).

For lattice fermions, multigrid targets critical slowing down caused by low-energy modes. In twisted-mass simulations, DD-$L^h u^h=f^h.$8AMG combines adaptive aggregation, domain-decomposition smoothing, and coarse-grid correction. It is used selectively with CG and multi-mass-shift CG: multigrid handles the most ill-conditioned systems, while CG or MMS-CG handles sufficiently well-conditioned or shifted systems. In physical-point Hybrid Monte Carlo simulations, the reported end-to-end speedups are approximately a factor of 8 for $L^h u^h=f^h.$9 and a factor of 5 for $L^h$0 (Bacchio et al., 2017).

Staggered fermions require a spectral transformation before Galerkin coarsening. Direct adaptive Galerkin multigrid applied to the staggered Dirac operator can generate spurious coarse low eigenvalues and amplify high-frequency errors through an oblique coarse-correction projector. The transformation

$L^h$1

where $L^h$2 is a block-local Kähler–Dirac operator, produces an approximately shifted-unitary spectrum. The transformed operator is then coarsened using adaptive near-null vectors, chiral doubling, and recursive K-cycles. In two-dimensional Schwinger-model tests, fine and intermediate iteration counts are approximately independent of mass and volume, while critical slowing down is transferred to the small coarsest problem (Brower et al., 2018).

Wilson-clover multigrid in Grid applies multigrid directly to the non-Hermitian Dirac operator and preserves $L^h$3-Hermiticity by chiral separation. Its coarse operator is

$L^h$4

The implementation supports MR, GMRES, FGMRES, recursive coarsening, mixed-precision preconditioning, SIMD vectorization, OpenMP, and MPI. An optimization of the restriction kernel produced approximately a 15-fold speedup for that kernel and approximately a twofold speedup of complete setup. The initial implementation had more expensive setup than DDalphaAMG but faster solves; subsequent improvements reduced total runtime to approximately 60% of DDalphaAMG on the tested lattice (Richtmann et al., 2019).

6. Data-driven and extended multigrid formulations

Several approaches reinterpret multigrid’s coarse representation or smoother without abandoning its error-elimination objective.

Mode Multigrid (MMG) applies dynamic mode decomposition to a sequence of flow solutions rather than constructing spatially coarsened meshes. DMD separates a dominant approximately zero-frequency steady mode from decaying oscillatory modes. Retaining the first mode and reconstructing the flow in physical space acts as a modal analogue of coarse-grid error filtering. The method leaves the physical mesh, discretization, and time integrator unchanged. Tests on unstructured finite-volume flows report iteration reductions of approximately 50%–80% and speedups of about $L^h$5–$L^h$6, including three-dimensional ONERA M6 and viscous anisotropic-mesh cases (Liu et al., 2018).

A learned prolongation can also replace conventional interpolation. A super-resolution GAN has been used as the coarse-to-fine operator in a two-dimensional pressure-Poisson multigrid solver. The GAN operates on $L^h$7 coarse tiles and produces $L^h$8 fine tiles. In the reported experiments, GAN-only interpolation improved mean convergence relative to cubic spline interpolation, while alternating GAN and spline interpolation converged faster for all 100 test fields. The method is distribution-dependent: it is trained on pressure fields from biperiodic incompressible-flow simulations, and no general stability or transfer-operator theory is established (Holguin et al., 2021).

Neural Multigrid Architectures formulate geometric multigrid as a matrix-free linear neural network. Convolutional and transposed-convolutional layers implement restriction, prolongation, residual computation, coarse correction, and smoothing. Polynomial smoothers are trainable because they require only repeated matrix-vector products. Training minimizes a stochastic estimate of the spectral radius of the error-propagation matrix,

$L^h$9

Parameter sharing and serialization enlarge the effective receptive field under grid refinement while keeping the number of trainable parameters independent of the number of levels. For tested elliptic systems, learned smoothers produce spectral radii approximately two to five times smaller than a Jacobi-smoothed linear multigrid baseline. Generalization is strongest when geometric transfer operators are fixed and smoothing is learned; jointly learning restriction and smoothing is less reliable (Fanaskov, 2024).

Multigrid can also operate on dimensions other than physical space. In PDE-constrained optimization under uncertainty, a collective smoother updates all uncertainty samples associated with a spatial node simultaneously. The local saddle-point system has dimension $(L^h)^T=L^h,\qquad x^TL^hx>0$0, where $(L^h)^T=L^h,\qquad x^TL^hx>0$1 is the number of uncertainty samples, but its solution costs $(L^h)^T=L^h,\qquad x^TL^hx>0$2. All nodal solves therefore cost

$(L^h)^T=L^h,\qquad x^TL^hx>0$3

which is linear in the full algebraic dimension. The spatial hierarchy is coarsened while the stochastic sample set remains fixed. Collective multigrid has been used as a stationary method, a GMRES preconditioner, and an inner solver in semismooth Newton and risk-averse Newton methods (Ciaramella et al., 2023).

A further extension addresses singular fractional diffusion. For two-sided fractional finite-volume-element discretizations on mapped graded meshes, the matrices are dense but Toeplitz-like. Their spectral structure is represented by

$(L^h)^T=L^h,\qquad x^TL^hx>0$4

A parameter-free GMRES preconditioner uses geometry-aware transfer and automatically estimated damped Jacobi relaxation. Power-graded meshes target boundary singularities using

$(L^h)^T=L^h,\qquad x^TL^hx>0$5

For symmetric and moderately nonsymmetric cases, reported iteration counts are nearly mesh-independent; robustness deteriorates in the strongly anisotropic, nearly skew-symmetric regime (Donatelli et al., 2022).

7. Complexity, limitations, and research directions

The ideal multigrid work model is $(L^h)^T=L^h,\qquad x^TL^hx>0$6 when each level has work proportional to its number of unknowns and the number of effective cycles is bounded. In practice, complexity depends on hierarchy depth, transfer sparsity, coarse-grid work, coefficient regularity, anisotropy, mesh quality, and communication. Structured local operators can achieve linear work per level, while Toeplitz-like fractional operators may require $(L^h)^T=L^h,\qquad x^TL^hx>0$7 matrix-vector products. Auxiliary and cluster-tree methods commonly report $(L^h)^T=L^h,\qquad x^TL^hx>0$8 construction or storage because of hierarchy depth.

The principal limitations recur across formulations:

• Coarse-grid quality: excessive coarsening can make smooth error unrepresentable.
• Smoother compatibility: pointwise relaxation is weak for anisotropy, strong coupling, or pathological local matrix entries.
• Transfer quality: poor interpolation can destroy approximation properties or amplify high-frequency error.
• Irregular geometry: unstructured, curved, nonmanifold, and unfitted domains complicate prolongation.
• Small coarse grids: accelerator occupancy and distributed parallelism deteriorate as levels become small.
• Coefficient variation: discontinuous or highly heterogeneous coefficients may require specialized coarsening or smoothing.
• Non-Hermitian structure: direct Galerkin projection can fail when spectral geometry and oblique projectors are unfavorable.
• Setup cost: adaptive bases, graph construction, parameterization, and coarse operators may dominate when only a few solves are required.
• Communication: strong scaling eventually becomes limited by ghost exchange, interlevel transfer, synchronization, and coarse-grid serialization.
• Learned components: data-driven interpolation and modal filtering depend on training distributions, snapshot selection, boundary treatment, and stability safeguards.
• Unproven generality: many reported convergence results are restricted to specific operators, geometries, discretizations, or simplified analytical models.

Current research directions include semi-structured and heterogeneous coarsening, non-Galerkin operators, improved high-order and surface transfer operators, GPU and multi-GPU implementations, compatible discretizations for divergence constraints, adaptive smoothers for low-quality cells and cut cells, stochastic coarsening, robust treatments of non-Hermitian operators, and learned but structure-preserving transfer operators. Across these developments, the central principle remains unchanged: relaxation must eliminate the error components it can efficiently resolve, while a suitably constructed coarse representation must remove the smooth components that otherwise cause slow convergence.