Methods for the 1D Parametric Eigenvalue Problem: A Comparison Framework

What is on this page. A catalogue of methods for solving the parametric eigenvalue problem (pEVP) in one parameter dimension, framed as a basis for later head-to-head numerical experiments against the adaptive-greedy algorithm. The source work (Alghamdi 2025) introduces the adaptive-greedy method but does not compare it against alternatives — this document fills that gap so that future numerical experiments can pin “the greedy algorithm is X% faster / more accurate than method Y under conditions Z” as measurable claims rather than assumptions. Companion pages on this site document the rebuild’s own measured properties; the open methodological questions are collected in a companion research note (not yet published).

Scope. 1D only (one parameter dimension, d = 1). Symmetric-definite generalised eigenvalue problem: K(µ) u = λ M(µ) u with K(µ), M(µ) ∈ ℝ^{n×n} sparse, symmetric, M(µ) positive definite, µ ∈ [µ_lo, µ_hi] ⊂ ℝ. Affine-in-coefficient assembly is the reference problems’ regime; some methods exploit that structure, others don’t.

1. The reference method: adaptive-greedy

This is the baseline every other method is compared against. It is not the first thing tried in the literature — it is a recent contribution (Alghamdi 2025) — but in this document it plays the role of the reference because that is the codebase the comparisons feed back into.

Idea. Sample µ adaptively on a sparse grid. At each edge of the grid, perform a posteriori certification: a matched eigenpair is certified if its eigenvector overlap with the partner at the other endpoint is essentially diagonal (off-diagonal projection mass below t_π, eigenvalue separation above t_λ). Edges that fail certification are bisected; the loop runs until no edge is flagged or a maximum refinement level is reached.

Cost model. Let C_eig be the cost of one cold eigensolve at the working n_dof. The greedy loop performs one eigensolve per node, one cost-matrix + Hungarian assignment + a posteriori matrix construction per edge. Eigensolves dominate; assignment and projection are O(k³) and O(k² · n_dof) respectively, both negligible for n_dof ≫ k.

Let N_nodes(L) be the number of nodes after L refinement levels. Total cost ≈ N_nodes(L) · C_eig plus a small linear-algebra overhead per edge.

Where it shines.

  • Crossings can be located with grid spacing only as fine as needed for the verifier to certify around them.
  • Smooth regions get sparse sampling automatically — no manual grid design.
  • The basis built along the way (the snapshot store) is reusable for post-hoc ROM evaluation.

Where it breaks / what to watch for.

  • The termination criterion is structural, not error-bounded: the loop stops when no edge is flagged, but no a-priori bound ties that fixed point to a target accuracy |λ_h − λ| ≤ ε. What termination actually certifies is examined in a companion research note (not yet published) on open methodological questions.
  • Window mode’s window-entry / window-exit labels are a rebuild-only interpretation layer, not a measurement — no curve identification exists behind them. With the maximum refinement level raised above its default, this can drive the loop to the cap on boundary geometry rather than on basis failure.
  • Tolerances (t_π, t_λ) are chosen by trial and error, with no principled selection rule; a 20% perturbation can flip termination level and degrade ROM accuracy by an order of magnitude — see the tolerance-perturbation page.
  • Reproducing the reference §3.3.3 Table 2D Error per-level counts failed across the full set of tunable knobs available to the rebuild (verdict mode, pruning rule, eigenvector cluster rewrite, cluster exemption, refinement topology). The published table appears to come from a driver not in the source dump; the reproduction gap is documented in a companion research note (not yet published).

2. Reference baselines

These are not comparison “methods” in the same sense as §3–§7. They are what greedy is measured against: the cost denominator (§2A) and the branch-identity oracle (§2B). They are the honest framing of “how cheap can we be without being clever” (§2A) and “what would a ground-truth tracker cost” (§2B).

2A. Uniform dense grid + cold eigensolve per node

Idea. Sample µ_j = µ_lo + j · h for j = 0..N. At each µ_j solve K(µ_j) v = λ M(µ_j) v from scratch.

Cost. N · C_eig exactly. No reuse.

Where it shines.

  • Trivial to implement and debug — this is the rebuild’s dense reference sweep used to generate truth data for accuracy measurements.
  • Embarrassingly parallel across µ.
  • Error from grid spacing is bounded by smoothness of λ_n(µ) and obvious to control: halve h, error roughly halves between crossings.
  • No tolerance trial-and-error; no termination ambiguity.

Where it breaks.

  • Sample count is fixed in advance — wasteful in smooth regions, potentially under-resolved near crossings.
  • Branch tracking is post-hoc (see 2B); the raw eigensolve output is sorted-ascending and does not respect curve identity.
  • No reusable basis emerges; downstream ROM evaluation requires a separate construction.

Role in comparison. §2A is the cost denominator for every method below — not a competitor of the greedy algorithm. Pick a target accuracy ε. Find the coarsest uniform grid such that max_µ |λ_uniform(µ) − λ_dense(µ)| < ε. The adaptivity ratio N_nodes_uniform(ε) / N_nodes_greedy(ε) on §3.3.2 is a number the reference leaves unreported, and the headline figure any greedy-vs- alternative experiment depends on. The question is what the ratio is, not which side wins.

2B. Dense grid + Hungarian branch tracking

Idea. As 2A, plus: between every adjacent pair (µ_j, µ_{j+1}), run cost-matrix + Hungarian assignment on eigenvector overlaps to tag each eigenpair with a continuous branch label. This is the rebuild’s branch-tracking pass used to assemble tracked-spectrum figures.

Cost. 2A plus O((N − 1) · k²) matching work — dominated by 2A’s eigensolves at any realistic problem size.

Where it shines. Solves the curve-identity problem that the greedy algorithm’s Window mode cannot resolve. Output is “λ_n(µ_j) for branch n = 1..k, j = 0..N” — a fully tracked spectrum.

Where it breaks. Cost is the cost of 2A; no reuse of parametric structure. Hungarian matching itself fails the same way Window-mode matching does when adjacent-sample eigenvectors collapse near a crossing — the symptom is gauge-dependent branch label swaps.

Role in comparison. §2B is the branch-identity oracle: the “expensive correct” baseline for branch tracking that the greedy is trying to be a cheaper version of. Useful both as ground-truth for any branch-identity metric and as a cost upper bound that any pEVP method should beat.

3. Solver-level continuation (per-eigensolve reuse)

These methods accelerate the individual eigensolve at each µ_j by reusing work from the previous one. Sample count stays the same as 2A; per-sample cost goes down.

3A. Warm-started LOBPCG along a µ-path

Idea. Use the eigenvectors at µ_{j-1} as the initial block for LOBPCG at µ_j. Optionally with Jacobi or AMG preconditioning.

Cost. Per-iteration ≈ a few matvecs; total per µ_j is O(k · n_dof) in matvecs plus O(k³) Rayleigh–Ritz, with iteration count dropping sharply as the warm start improves.

Where it shines. Smooth regions where consecutive eigenpairs are close — convergence in 2-5 iterations.

Where it breaks.

  • In the rebuild’s measurements at n_dof ≤ 10⁴ (the reference problems’ regime), warm-started LOBPCG with Jacobi preconditioning is 5–10× slower than ARPACK shift-invert. The crossover would require either much larger n_dof or AMG-class preconditioning.
  • Block reorganization at crossings is awkward — the closest-eigenvalue continuation isn’t always the right pairing.

Head-to-head against the greedy. The experiment runs 2A with warm-started LOBPCG at each node and compares total wall time and accuracy against greedy’s N_nodes · C_eig. It measures whether solver-level reuse alone can close the sample-count gap that adaptive sampling exploits. Caveat: needs an n_dof large enough that LOBPCG is competitive in the first place — the reference problems’ regime isn’t.

3B. Krylov–Schur with subspace recycling

Idea. Same family as 3A; SLEPc-style implementation. Reuse the prior step’s Krylov subspace as the starting basis.

Cost. Comparable to 3A; the trade-off between block size and restart depth tunes differently.

Status. Not implemented in the rebuild. Worth measuring once AMG-class preconditioning is integrated; until then ARPACK cold-start is the default for principled reasons.

4. Smooth-regime perturbation / ODE continuation

These methods don’t sample the parameter at all in smooth regions — they integrate the eigenpair forward in µ using analytic derivatives.

4A. Analytic eigenvalue derivative (Rellich / Hellmann–Feynman)

Idea. For a simple eigenpair (λ(µ), v(µ)) of a symmetric problem:

dλ/dµ = vᵀ (∂K/∂µ − λ ∂M/∂µ) v / (vᵀ M v)

The eigenvector derivative satisfies a singular linear system that can be solved by either spectral projection or by augmenting with the normalisation constraint. For affine-in-µ problems (the reference problems’ regime), ∂K/∂µ and ∂M/∂µ are constant — assembled once.

Cost per step. One Taylor (or Padé) extrapolation; no eigensolve required. Per-step cost is one matvec for the derivative formula plus the linear solve for v'. Order of magnitude cheaper than a full eigensolve.

Where it shines. Smooth, well-separated branches. The predictor’s accuracy is O(h^p) for a p-th-order Taylor step; stepping all the way across [µ_lo, µ_hi] from one cold solve at µ_lo is feasible when no crossings occur.

Where it breaks.

  • Essential singularity at every crossing. dλ_i/dµ blows up as λ_i → λ_j because the eigenvector derivative requires (K − λ M) to be invertible on the orthogonal complement of v, and that fails at crossings. Two distinct branches need to exchange eigenvectors smoothly, which the formula above does not describe.
  • Branch identity is fragile: which eigenpair (λ_i, v_i) the integrator “is on” is determined by the initial pick; near an avoided crossing the integrator may follow the wrong branch.
  • Requires a corrector (Newton or a full eigensolve) periodically to prevent drift — pure predictor accumulates error.

Head-to-head against the greedy. Measures the cost-per-unit-µ of analytic continuation in the smooth-only regime, and how far the method degrades when crossings are present. The experiment benchmarks 4A on a problem with no crossings (or with crossings manually excised) and on the §3.3.2 problem with crossings; the crossing case answers whether a pure predictor without branch switching is salvageable at all.

4B. Predictor–corrector continuation (arc-length)

Idea. 4A as predictor; full eigensolve at fixed µ (or arc-length-parameterised step) as corrector. Branch switching at detected crossings via a test function (e.g., sign change of an eigenvalue gap).

Cost. Per accepted step: one predictor (cheap) + one corrector (expensive but possibly warm-started). Step length adapts to local curvature of λ(µ).

Where it shines. Smooth-regime efficiency of 4A combined with recovery at crossings. This is the classical numerical continuation playbook (Keller, Allgower–Georg).

Where it breaks.

  • Branch switching at crossings is delicate. The greedy algorithm’s Hungarian matching is, in effect, a discrete branch-switching rule applied at every sampled edge. 4B has to do the same thing but with finer control over when to apply it (only at detected crossings, not at every step).
  • Corrector cost dominates when crossings are dense.

Head-to-head against the greedy. Measures where each method spends its eigensolves: 4B does eigensolves at detected crossings; the greedy does eigensolves at every bisection level around crossings, including levels that fire on window-entry / window-exit interpretation rather than on basis failure. The metric is eigensolve count to a fixed ε accuracy on §3.3.2.

4C. Subspace continuation on the Grassmannian

Idea. Instead of tracking individual eigenpairs, track the invariant subspace spanned by the lowest k eigenvectors. The subspace is robust through crossings — only the basis of the subspace rotates. Methods: Mehrmann–Schwetlick, Davis–Kahan perturbation bounds.

Cost. A k-dimensional subspace ODE integrated in µ; periodic projection back onto the eigenproblem manifold via Rayleigh–Ritz.

Where it shines. Robust through crossings — the bottleneck of 4A and 4B. The price is that individual branches are not tracked during the continuation; they are recovered post-hoc by diagonalising the projected matrix.

Where it breaks. Mathematically clean, computationally non-trivial; not a “turn the crank” method. The literature is specialised and the off-the-shelf software is thin.

Head-to-head. The most interesting comparison for the crossing-heavy regime the adaptive-greedy method targets. Worth a literature survey before any numerical experiment.

5. Surrogate / reduced-basis methods (offline–online split)

These methods build a global surrogate from a training set of snapshots, then evaluate cheaply at any query µ.

5A. Polynomial / Chebyshev interpolation of K(µ), M(µ)

Idea. Sample K(µ_j), M(µ_j) at Chebyshev nodes; interpolate component-wise. The resulting K̃(µ), M̃(µ) are polynomial in µ. The reduced problem K̃(µ) u = λ M̃(µ) u is a polynomial eigenvalue problem (PEP) — solvable once globally via companion linearisation.

Cost. N_cheb assemblies offline; one PEP solve online at the chosen µ (no extra eigensolves for additional µ).

Where it shines. Affine-in-µ problems: interpolation is exact with N_cheb = degree(µ) + 1 nodes. The §3.3.2 benchmark problem has K(µ) = µ⁻² K_1 + K_2 — not polynomial in µ, but polynomial in µ⁻². So Chebyshev in t = µ⁻² gives an exact representation with 2 nodes.

Where it breaks. λ(µ) is not polynomial in µ even when K(µ), M(µ) are — it’s an algebraic function. The PEP-of-an- interpolant approach gives access to all n_dof eigenvalues, not just the lowest k, so band-selection (first-k or window) must be imposed afterwards.

Head-to-head against the greedy. For affine problems, this measures whether the greedy’s sample count exceeds the structural minimum dictated by the polynomial degree of K(µ), M(µ) (2 for §3.3.2). The reference L=0..L=3 grid sizes for §3.3.2 are far above that minimum; the experiment quantifies how far above and allocates the excess between branch-identity costs and verdict-loop overhead.

5B. Reduced-basis (POD on snapshots)

Idea. Sample µ at a training set; collect eigenvectors; SVD → reduced basis V ∈ ℝ^{n_dof × r}. For any query µ, solve the reduced problem Vᵀ K(µ) V x = λ Vᵀ M(µ) V x.

Cost. Offline: N_train · C_eig + SVD. Online: one r × r dense GEVP per query.

Where it shines. Many queries on a fixed parameter range. Same idea as the snapshot store in the rebuild — but with fixed offline training, not adaptive.

Where it breaks.

  • Choice of training set is manual or heuristic. The greedy algorithm exists precisely to remove this choice.
  • Accuracy depends on whether the training set captured the parametric behaviour; no a posteriori certificate.

Head-to-head against the greedy. Train 5B on the greedy’s terminal grid versus on a uniform grid of the same size, and compare ROM accuracy at unseen µ. This measures whether the greedy’s sample selection adds value beyond just selecting that many samples. If a uniform grid of the same node count gives comparable accuracy, the greedy’s adaptivity is buying very little.

5C. Greedy reduced basis (Maday–Patera RBM)

Idea. Classic reduced-basis method with a greedy training-set selection driven by an a posteriori error estimator. The adaptive-greedy algorithm sits in this family — the specific innovations being (i) sparse-grid topology over parameter space, (ii) the projection-mass verdict, and (iii) window-entry / window-exit labelling for windowed band selection.

Head-to-head against the rebuild’s specific algorithm. Substitute the reference method’s verdict for the classical residual-norm estimator η(µ) = ‖K(µ) V x − λ M(µ) V x‖ / ‖M(µ) V x‖ and measure how the refinement trajectory differs. The classical estimator is error-driven; the adaptive-greedy’s is projection-driven. These are not equivalent, and the experiment quantifies which one selects more useful samples on a fixed problem.

5D. Direct interpolation of λ_n(µ) per branch

Idea. Sample a dense µ-grid, run 2B to get branch-tracked eigenvalues, then fit a smooth interpolant (cubic spline, Chebyshev, radial basis function) to each branch independently.

Cost. Same as 2B offline; trivial online.

Where it shines. Smooth branches in smooth regions; the interpolant captures λ_n(µ) to machine precision with few nodes per branch.

Where it breaks. Crossings: interpolating across a crossing mixes branches. Has to be combined with explicit crossing detection to subdivide the interval.

Head-to-head. This is the simplest possible surrogate that uses branch information. If the goal is just λ_n(µ) evaluation (not eigenvector evaluation), the experiment measures whether anything fancier than per-branch interpolation pays off in 1D.

6. Contour-integral spectral methods

These methods are window-mode native: they compute eigenpairs in a prescribed [λ_lo, λ_hi] without computing the rest of the spectrum.

6A. FEAST / Sakurai–Sugiura at each µ_j

Idea. Enclose [λ_lo, λ_hi] in a complex contour. Approximate the spectral projector onto the enclosed eigenspace by Gauss–Legendre quadrature of the resolvent (zM − K)⁻¹M over the contour. Diagonalise the projector to recover the enclosed eigenpairs.

Cost. Per µ_j: one linear solve at each quadrature node (≈ 8–16 nodes typical). For affine K(µ), the per-node factorisations can be precomputed once across all µ if the resolvent’s µ-dependence is handled symbolically. Otherwise per-µ cost is O(N_quad · C_solve).

Where it shines.

  • Window mode native — solves exactly the band-selection problem the adaptive-greedy’s Window mode tries to handle (whose window-entry / window-exit semantics remain ambiguous in the greedy).
  • Embarrassingly parallel across both µ and quadrature nodes.
  • No “first-k” cutoff issue; no truncation-gauge problem.

Where it breaks.

  • Number of eigenpairs in the window is not known in advance; FEAST needs a guess for the subspace dimension. Too small → missed eigenvalues; too large → wasted work.
  • Conditioning of the contour quadrature depends on the gap between the contour and the nearest eigenvalues outside the window.

Head-to-head against greedy Window mode. Solve §3.3.2’s window problem ([0, 270] per the reference choice) by FEAST at each µ on a uniform grid, and compare against the greedy algorithm in Window mode. The experiment measures the trade-off between two qualitatively different cost axes: FEAST has no window-entry / window-exit ambiguity by construction but pays the full per-µ cost; the greedy reduces sample count but inherits that interpretive ambiguity at termination.

7. Symbolic / characteristic-polynomial methods (very low-dim only)

For very small n_dof (say, n_dof ≤ 20), the characteristic polynomial p(µ, λ) = det(K(µ) − λ M(µ)) can be computed symbolically in µ if K(µ), M(µ) are polynomial (or rational) in µ. The zero set of p(µ, λ) is an algebraic curve in (µ, λ)-space, and the branches are the connected components of that curve.

Cost. Symbolic determinant: O(n³) in symbolic operations, exponential in the polynomial degree of K(µ), M(µ). Solving p(µ_j, λ) = 0 at a fixed µ_j is a univariate root-finding in λ.

Status. Not practical at the reference problems’ n_dof ≈ 10⁴. Worth mentioning because it is the mathematically exact answer in the toy regime where it applies — useful as a sanity check on implementations at, e.g., n_dof = 4.

Head-to-head. A unit-test-scale verification, not a benchmark.

8. Comparison axes — what an experiment would measure

Any head-to-head experiment should report a fixed tuple, computed against the same dense reference (a uniform fine-mesh sweep, §2A):

Axis Definition Why it matters
Eigensolve count Number of eigensolver calls The dominant cost at any realistic n_dof. Sample count × per-sample solve count.
Wall time Total compute time to reach a fixed accuracy ε Captures secondary costs (matching, projection, assembly) that the eigensolve count ignores.
λ-accuracy max_µ max_n |λ_method(µ, n) − λ_truth(µ, n)| after Hungarian-aligning branches The actual quantity the user cares about.
Eigenvector accuracy max_µ max_n subspace angle between cluster-equivalent eigenspaces The quantity downstream physics computations care about.
Branch-identity error rate Fraction of (µ, n) pairs where the method’s branch label disagrees with the dense reference’s Measures whether the method actually solves “track curves” or only “find eigenvalues.” Greedy Window mode answers this poorly — no curve identification exists behind its window-entry / window-exit labels.
Crossing detection True / false positives / negatives on crossings, vs. dense reference Specific to crossing-handling claims. The greedy’s refine-or-certify verdict is silent on crossings; window-entry / window-exit labels claim them but cannot certify them.
Robustness to tolerance perturbation Same metrics under ±20% perturbation of all tolerances / step sizes The reference method’s tolerances are trial-and-error, with no principled selection rule; robustness is a first-class property, not an afterthought.

Each experiment fixes a method, a target accuracy ε, and reports the seven-tuple. The reference method (greedy, §1) provides the baseline column; each other method provides one or more columns.

9. Candidate experiments worth scoping

Three experiments stand out as the most informative if the comparison framework is ever exercised end-to-end. None has been committed to; this list scopes the questions, not the build.

  1. Adaptivity ratio on §3.3.2. Greedy vs. uniform dense grid (§2A) on the reference 1D problem. Question: at the same ε, what is N_nodes_uniform(ε) / N_nodes_greedy(ε)? This is the headline number the reference leaves unreported, and the denominator every other comparison in this catalogue feeds into.

  2. Greedy vs. predictor–corrector with branch switching (§4B). Same problem. Question: does locating crossings explicitly (test-function detection + branch switch) outperform locating them implicitly (bisection until certification fires)? This is the most direct comparison of the adaptive-greedy method’s contribution against the classical numerical-continuation playbook.

  3. Greedy Window mode vs. FEAST per-µ on a uniform grid (§6A). Specifically on §3.3.2 with the reference I_λ = [0, 270]. Question: does the greedy’s sample-count advantage outweigh FEAST’s semantic-clarity advantage (no window-entry / window-exit ambiguity, no truncation gauge)?

Experiment #1 is the most consequential — it answers whether the greedy is worth its complexity over a dense grid at all, and calibrates the denominator the other two are reported against.

The continuation methods catalogued in §§3–4 are also the candidate attacks on the 6D wall-time wall — after the snapshot store ceased to be the binding constraint, the §4.4.5 6D run is capped by the sheer number of full eigensolves, not by memory. That feasibility problem — is there a crossing-safe per-step continuation cheap enough to replace the cold-start eigensolve at each new parameter point and finish 6D in practical wall time? — is posed in a companion research note (not yet published) on open problems, which draws on this section’s cost / crossing-safety trade-offs.

10. What this document is not

Not a literature survey (each method gets one paragraph, not its canonical references), not a benchmark (no numbers — the point is to fix the axes before running anything), and not exhaustive (block iterative methods, trace-based spectral indicators, random-projection methods, and FE multigrid eigensolvers are explicitly omitted).

Last updated: 2026-06-08.