GD–λAMG: Multigrid-Accelerated Davidson

• GD–λAMG is an eigensolver that combines the generalized Davidson method with a multigrid-accelerated inner iteration to compute extremal eigenpairs of large, sparse Hermitian matrices in lattice QCD.
• It features dynamic updates of coarse operator components and leverages local coherence to maintain preconditioning efficiency while solving correction equations.
• Benchmark studies indicate that GD–λAMG achieves significant speedups—up to 10× faster than traditional methods—ensuring accurate spectral analysis for improved QCD simulations.

Multigrid-accelerated Davidson (GD–λAMG) refers to a class of eigensolvers that combine the generalized Davidson method with a multigrid-accelerated inner iteration for efficiently computing extremal (typically smallest magnitude) eigenvalues and eigenvectors of large, sparse Hermitian matrices. This approach is particularly relevant in the spectral analysis of the Hermitian Wilson–Dirac operator in lattice quantum chromodynamics (QCD), where resolving the low-lying spectrum is crucial for various physical observables and for accelerating linear system solutions via deflation (Frommer et al., 2020).

1. Problem Formulation and Context in Lattice QCD

The Hermitian Wilson–Dirac operator $Q$ on a four-dimensional lattice $\mathcal{L}$ of size $N_t\times N_s^3$ with quark mass $m_0$ is constructed as

$$D\psi(x) = (m_0 + 4/a) \psi(x) - \frac{1}{2a} \sum_{\mu=0}^3 [ (I - \gamma_\mu)\otimes U_\mu(x)]\psi(x+\hat\mu) - \frac{1}{2a} \sum_{\mu=0}^3 [ (I + \gamma_\mu)\otimes U^\dagger_\mu(x-\hat\mu)] \psi(x-\hat\mu)$$

with $\Gamma_5 = \gamma_0\gamma_1\gamma_2\gamma_3$. The Hermitian operator is defined as $Q := \Gamma_5 D$, so $Qx = \lambda x$ becomes the target eigenproblem with $x\in\mathbb{C}^N$, $\lambda\in\mathbb{R}$, specifically focusing on eigenpairs with small $\mathcal{L}$0.

Small$\mathcal{L}$1 eigenpairs of $\mathcal{L}$2 are essential in noise-reduction, all-to-all propagators, low-mode averaging, linear solve deflation, and in the study of spontaneous chiral symmetry breaking via the Banks–Casher relation. The spectral density near zero directly informs physical phenomena and the efficiency of many computational techniques (Frommer et al., 2020).

2. Generalized Davidson Method: Core Algorithmic Components

The generalized Davidson method forms the outer iteration framework. An orthonormal basis $\mathcal{L}$3 is maintained, spanning the current search space. At each step, harmonic Ritz pairs $\mathcal{L}$4 are found by the generalized eigenproblem

$\mathcal{L}$5

producing approximate eigenvalues $\mathcal{L}$6 and corresponding Ritz vectors $\mathcal{L}$7.

The residual $\mathcal{L}$8 quantifies the approximation error. A new search direction $\mathcal{L}$9 is obtained by approximately solving the shifted equation (the correction equation)

$N_t\times N_s^3$0

with the shift $N_t\times N_s^3$1 typically chosen close to $N_t\times N_s^3$2. The new vector $N_t\times N_s^3$3 is orthogonalized against the current basis and any locked (converged) eigenvectors to expand the search space. This process is repeated until $N_t\times N_s^3$4 target eigenpairs are identified.

3. Multigrid-Accelerated DD–$N_t\times N_s^3$5AMG Inner Iterations

At each Davidson step, the correction equation takes the form $N_t\times N_s^3$6 with shift $N_t\times N_s^3$7, where $N_t\times N_s^3$8 is the residual. The algorithm employs left preconditioning with $N_t\times N_s^3$9, yielding $m_0$0, with $m_0$1.

A multigrid hierarchy is constructed using aggregation-based interpolation $m_0$2, built from $m_0$3 “test vectors” (approximations to the smallest$m_0$4 modes). The restriction is $m_0$5, and the coarse operator is $m_0$6. For the shifted system, $m_0$7 with $m_0$8 and $m_0$9 is used.

A two-level error propagator for $$D\psi(x) = (m_0 + 4/a) \psi(x) - \frac{1}{2a} \sum_{\mu=0}^3 [ (I - \gamma_\mu)\otimes U_\mu(x)]\psi(x+\hat\mu) - \frac{1}{2a} \sum_{\mu=0}^3 [ (I + \gamma_\mu)\otimes U^\dagger_\mu(x-\hat\mu)] \psi(x-\hat\mu)$$0 is constructed as

$$D\psi(x) = (m_0 + 4/a) \psi(x) - \frac{1}{2a} \sum_{\mu=0}^3 [ (I - \gamma_\mu)\otimes U_\mu(x)]\psi(x+\hat\mu) - \frac{1}{2a} \sum_{\mu=0}^3 [ (I + \gamma_\mu)\otimes U^\dagger_\mu(x-\hat\mu)] \psi(x-\hat\mu)$$1

with $$D\psi(x) = (m_0 + 4/a) \psi(x) - \frac{1}{2a} \sum_{\mu=0}^3 [ (I - \gamma_\mu)\otimes U_\mu(x)]\psi(x+\hat\mu) - \frac{1}{2a} \sum_{\mu=0}^3 [ (I + \gamma_\mu)\otimes U^\dagger_\mu(x-\hat\mu)] \psi(x-\hat\mu)$$2 a smoother (e.g., SAP or GMRES) and $$D\psi(x) = (m_0 + 4/a) \psi(x) - \frac{1}{2a} \sum_{\mu=0}^3 [ (I - \gamma_\mu)\otimes U_\mu(x)]\psi(x+\hat\mu) - \frac{1}{2a} \sum_{\mu=0}^3 [ (I + \gamma_\mu)\otimes U^\dagger_\mu(x-\hat\mu)] \psi(x-\hat\mu)$$3 specifying pre- and post-smoothing steps. The approximate inverse $$D\psi(x) = (m_0 + 4/a) \psi(x) - \frac{1}{2a} \sum_{\mu=0}^3 [ (I - \gamma_\mu)\otimes U_\mu(x)]\psi(x+\hat\mu) - \frac{1}{2a} \sum_{\mu=0}^3 [ (I + \gamma_\mu)\otimes U^\dagger_\mu(x-\hat\mu)] \psi(x-\hat\mu)$$4 is realized by a single V-cycle: pre-smoothing, coarse correction, post-smoothing.

4. Dynamic Coarse Operator Updates and Local Coherence

“Local coherence” refers to the property that small$$D\psi(x) = (m_0 + 4/a) \psi(x) - \frac{1}{2a} \sum_{\mu=0}^3 [ (I - \gamma_\mu)\otimes U_\mu(x)]\psi(x+\hat\mu) - \frac{1}{2a} \sum_{\mu=0}^3 [ (I + \gamma_\mu)\otimes U^\dagger_\mu(x-\hat\mu)] \psi(x-\hat\mu)$$5 eigenvectors exhibit similar structure on each aggregate. For a vector $$D\psi(x) = (m_0 + 4/a) \psi(x) - \frac{1}{2a} \sum_{\mu=0}^3 [ (I - \gamma_\mu)\otimes U_\mu(x)]\psi(x+\hat\mu) - \frac{1}{2a} \sum_{\mu=0}^3 [ (I + \gamma_\mu)\otimes U^\dagger_\mu(x-\hat\mu)] \psi(x-\hat\mu)$$6, the coherence metric $$D\psi(x) = (m_0 + 4/a) \psi(x) - \frac{1}{2a} \sum_{\mu=0}^3 [ (I - \gamma_\mu)\otimes U_\mu(x)]\psi(x+\hat\mu) - \frac{1}{2a} \sum_{\mu=0}^3 [ (I + \gamma_\mu)\otimes U^\dagger_\mu(x-\hat\mu)] \psi(x-\hat\mu)$$7 is used; values near 1 indicate that the interpolation $$D\psi(x) = (m_0 + 4/a) \psi(x) - \frac{1}{2a} \sum_{\mu=0}^3 [ (I - \gamma_\mu)\otimes U_\mu(x)]\psi(x+\hat\mu) - \frac{1}{2a} \sum_{\mu=0}^3 [ (I + \gamma_\mu)\otimes U^\dagger_\mu(x-\hat\mu)] \psi(x-\hat\mu)$$8 provides an accurate representation.

If the multigrid interpolation $$D\psi(x) = (m_0 + 4/a) \psi(x) - \frac{1}{2a} \sum_{\mu=0}^3 [ (I - \gamma_\mu)\otimes U_\mu(x)]\psi(x+\hat\mu) - \frac{1}{2a} \sum_{\mu=0}^3 [ (I + \gamma_\mu)\otimes U^\dagger_\mu(x-\hat\mu)] \psi(x-\hat\mu)$$9 is constructed from the smallest$\Gamma_5 = \gamma_0\gamma_1\gamma_2\gamma_3$0 eigenmodes, as the Davidson method advances and focuses on eigenvalues further from zero, $\Gamma_5 = \gamma_0\gamma_1\gamma_2\gamma_3$1 drops, resulting in reduced efficiency of the inner multigrid solve. GD–$\Gamma_5 = \gamma_0\gamma_1\gamma_2\gamma_3$2AMG addresses this by dynamically updating $\Gamma_5 = \gamma_0\gamma_1\gamma_2\gamma_3$3 and $\Gamma_5 = \gamma_0\gamma_1\gamma_2\gamma_3$4. When new eigenvectors $\Gamma_5 = \gamma_0\gamma_1\gamma_2\gamma_3$5 (corresponding to converged eigenvalues $\Gamma_5 = \gamma_0\gamma_1\gamma_2\gamma_3$6) become available, the $\Gamma_5 = \gamma_0\gamma_1\gamma_2\gamma_3$7 vectors with eigenvalues closest to the current target $\Gamma_5 = \gamma_0\gamma_1\gamma_2\gamma_3$8 are selected to rebuild $\Gamma_5 = \gamma_0\gamma_1\gamma_2\gamma_3$9—specifically,

$Q := \Gamma_5 D$0

This ensures that local coherence remains high throughout the Davidson iterations, maintaining the efficacy of the multigrid acceleration (Frommer et al., 2020).

5. Detailed Algorithmic Workflow (GD–λAMG)

The overall algorithm proceeds as follows:

1. Initialize $Q := \Gamma_5 D$1 with random/test vectors; build $Q := \Gamma_5 D$2, $Q := \Gamma_5 D$3; set $Q := \Gamma_5 D$4, $Q := \Gamma_5 D$5.
2. For each $Q := \Gamma_5 D$6:
   • Orthogonalize new guess $Q := \Gamma_5 D$7 to $Q := \Gamma_5 D$8 and $Q := \Gamma_5 D$9; normalize to obtain $Qx = \lambda x$0; augment $Qx = \lambda x$1.
   • Compute harmonic Ritz pairs $Qx = \lambda x$2 from $Qx = \lambda x$3; pick $Qx = \lambda x$4 closest to 0 not yet locked.
   • Set $Qx = \lambda x$5, $Qx = \lambda x$6.
   • If $Qx = \lambda x$7, lock $Qx = \lambda x$8: $Qx = \lambda x$9, $x\in\mathbb{C}^N$0; select $x\in\mathbb{C}^N$1 locked vectors nearest to next $x\in\mathbb{C}^N$2; rebuild $x\in\mathbb{C}^N$3, $x\in\mathbb{C}^N$4; restart subspace (size $x\in\mathbb{C}^N$5).
   • Solve preconditioned correction by one V-cycle of DD–$x\in\mathbb{C}^N$6AMG: $x\in\mathbb{C}^N$7.
3. Repeat until $x\in\mathbb{C}^N$8 eigenpairs are obtained.

This process tightly couples the outer eigensolver with the inner multigrid-accelerated correction step, with frequent adaptation of the coarse-level components.

6. Performance Characteristics and Comparative Benchmarks

Performance tests on large lattices, e.g., $x\in\mathbb{C}^N$9, show that GD–$\lambda\in\mathbb{R}$0AMG achieves substantial speedups over conventional eigensolvers. Key findings include:

• $\lambda\in\mathbb{R}$1-preconditioning yields an approximate $\lambda\in\mathbb{R}$2 reduction in solve time (e.g., 83 core-h to 41 core-h for 100 eigenpairs).
• DD–$\lambda\in\mathbb{R}$3AMG preconditioning reduces inner iteration counts per eigenvalue from $\lambda\in\mathbb{R}$4–$\lambda\in\mathbb{R}$5 (for pure GMRESR) to $\lambda\in\mathbb{R}$6 (for one V-cycle), corresponding to a $\lambda\in\mathbb{R}$7 speed-up.
• Without interpolation/coarse-grid updates, cost scales as $\lambda\in\mathbb{R}$8 for $\lambda\in\mathbb{R}$9 eigenpairs; with updates, cost is near-linear in $\mathcal{L}$00, halving total time at $\mathcal{L}$01.
• Benchmarking against PARPACK (with optimal polynomial filter) and PRIMME (with the same AMG), GD–$\mathcal{L}$02AMG is up to $\mathcal{L}$03 faster than PARPACK and $\mathcal{L}$04–$\mathcal{L}$05 faster than PRIMME (speedup increases with lattice size and $\mathcal{L}$06).
• For $\mathcal{L}$07 between $\mathcal{L}$08 and $\mathcal{L}$09, GD–$\mathcal{L}$10AMG maintains near-linear scaling, with the performance advantage over PRIMME increasing from $\mathcal{L}$11 to over $\mathcal{L}$12 (Frommer et al., 2020).

7. Spectral Analysis and Physical Relevance

Computed low-lying eigenvalues enable construction of the spectral density $\mathcal{L}$13, a key observable in lattice QCD. For a $\mathcal{L}$14 lattice with lattice spacing $\mathcal{L}$15 and $\mathcal{L}$16, the density from the smallest 1000 eigenvalues of $\mathcal{L}$17 is compared to the Banks–Casher relation, which in the continuum (for $\mathcal{L}$18, where $\mathcal{L}$19 is the bare quark mass, $\mathcal{L}$20) predicts

$\mathcal{L}$21

with $\mathcal{L}$22 the chiral condensate. The numerical results display nonzero density below $\mathcal{L}$23 and a smooth interpolation near $\mathcal{L}$24 rather than the divergence expected in the infinite-volume, continuum limit. This quantifies the impact of lattice discretization and volume artifacts in practical QCD simulations (Frommer et al., 2020).

In summary, GD–$\mathcal{L}$25AMG couples harmonic-Davidson iteration with dynamically updated DD–$\mathcal{L}$26AMG inner solves, exploiting local coherence and preconditioning to yield up to order-of-magnitude acceleration and near-linear scaling in both number of eigenpairs and lattice size, while preserving physically accurate spectral information critical for modern lattice QCD.