Generalized Eigenvalue Stabilization (GEVS)

• GEVS is a collection of techniques that stabilize generalized eigenvalue problems by controlling the spectrum and regularizing matrices such as B.
• It employs methods like regularization, whitening, and shift selection to suppress nonphysical eigenmodes and enhance numerical convergence.
• GEVS finds applications in explicit dynamics, quantum algorithms, and optimization, ensuring meaningful spectra and improved computational stability.

Generalized Eigenvalue Stabilization (GEVS) denotes a family of techniques used to make generalized eigenvalue problems, and computations driven by them, numerically stable, well posed, and spectrally interpretable. In its most classical form, the underlying problem is the generalized eigenvalue equation

$A\phi_i=\lambda_i B\phi_i,$

or, in matrix form, $A\Phi=B\Phi\Lambda$, with stabilization becoming necessary when $B$ is ill-conditioned, singular, or semidefinite, or when discretization artifacts generate nonphysical eigenmodes (Ghojogh et al., 2019). In more application-specific settings, GEVS also refers to selective spectral modifications of local pencils such as $K\varphi=\lambda M\varphi$ in explicit dynamics, where the objective is to suppress spectral outliers and recover a practical critical time step size (Bürchner et al., 9 Sep 2025). Across the literature, the common theme is explicit control of the spectrum, residuals, or invariant subspaces so that computed eigenpairs correspond to meaningful nearby problems, or so that problematic modes are neutralized without destroying the useful part of the spectrum (Stewart, 2024).

1. Mathematical setting and problem classes

The generalized eigenvalue problem (GEP) for symmetric matrices $A,B\in\mathbb{R}^{d\times d}$ is defined by

$$A\phi_i=\lambda_i B\phi_i,\qquad i\in\{1,\dots,d\},$$

with the standard eigenvalue problem recovered when $B=I$. The same structure arises as the first-order optimality condition of constrained quadratic optimization, both in vector form,

$$\underset{\phi}{\text{maximize}}\ \phi^\top A\phi \quad\text{subject to}\quad \phi^\top B\phi=1,$$

and in matrix form,

$$\underset{\Phi}{\text{maximize}}\ \mathrm{tr}(\Phi^\top A\Phi) \quad\text{subject to}\quad \Phi^\top B\Phi=I.$$

Equivalently, the generalized Rayleigh-Ritz quotient

$R(A,B;x)=\frac{x^\top Ax}{x^\top Bx}$

leads to the same stationarity condition (Ghojogh et al., 2019).

This algebraic structure recurs in disparate domains. In Fisher discriminant analysis one obtains $A\Phi=B\Phi\Lambda$0, while in kernel supervised principal component analysis the solution also requires a GEP (Ghojogh et al., 2019). In lattice gauge theory, the correlation-matrix formulation

$A\Phi=B\Phi\Lambda$1

is used to extract energies and matrix elements, with $A\Phi=B\Phi\Lambda$2 at large $A\Phi=B\Phi\Lambda$3 (0808.1017). In immersed explicit dynamics, the physically relevant spectrum is that of

$A\Phi=B\Phi\Lambda$4

where cut elements can create spurious high eigenvalues that collapse the critical time step (Bürchner et al., 9 Sep 2025).

This range of formulations suggests that GEVS is not a single algorithm but a spectral-control paradigm. The target may be invertibility of $A\Phi=B\Phi\Lambda$5, robustness of iterative eigensolvers, certification of computed eigenvalues, or suppression of discretization-induced outliers, but the organizing object remains a generalized eigenstructure.

2. Regularization, whitening, and preconditioning

A basic stabilization route is to regularize the constraint matrix $A\Phi=B\Phi\Lambda$6. If $A\Phi=B\Phi\Lambda$7 is invertible, the “quick and dirty” reduction

$A\Phi=B\Phi\Lambda$8

converts the problem to a standard eigenvalue problem. If $A\Phi=B\Phi\Lambda$9 is singular, the tutorial literature recommends adding a regularization term,

$B$0

with $B$1, so that $B$2 becomes invertible. In that setting, regularization both stabilizes the solution and makes the problem well posed (Ghojogh et al., 2019).

A more rigorous route is whitening. One first diagonalizes $B$3,

$B$4

then defines

$B$5

regularizing $B$6 by $B$7 when zero or very small eigenvalues are present. The transformed matrix

$B$8

is then solved as a standard eigenproblem, and the original eigenvectors are recovered via $B$9 (Ghojogh et al., 2019). The tutorial explicitly presents this whitening-based procedure as numerically robust and as a foundation for GEVS methods.

Recent large-scale optimization work reframes stabilization at the objective level. Instead of maximizing the classical Rayleigh quotient, which is described as numerically unstable in large-scale or ill-conditioned settings, a difference-based formulation minimizes

$K\varphi=\lambda M\varphi$0

Its Hessian satisfies a bounded positive curvature condition, $K\varphi=\lambda M\varphi$1 with $K\varphi=\lambda M\varphi$2, which yields global convergence with fixed stepsizes and no line search. In the same work, transform-domain preconditioning is based on a factor $K\varphi=\lambda M\varphi$3 such that

$K\varphi=\lambda M\varphi$4

leading to the preconditioned mirror descent update

$K\varphi=\lambda M\varphi$5

With $K\varphi=\lambda M\varphi$6 chosen as a Cholesky factor and $K\varphi=\lambda M\varphi$7, this update is algebraically equivalent to the power method (Liu et al., 3 Jul 2025).

A recurrent misconception is that stabilization is exhausted by replacing the GEP with $K\varphi=\lambda M\varphi$8. The literature instead treats the conditioning of $K\varphi=\lambda M\varphi$9, the geometry of the transformed problem, and the bias-variance tradeoff induced by $A,B\in\mathbb{R}^{d\times d}$0 as central design variables (Ghojogh et al., 2019).

3. Shift selection, perturbation, and spectral certification

For dense symmetric semidefinite GEPs with $A,B\in\mathbb{R}^{d\times d}$1 and $A,B\in\mathbb{R}^{d\times d}$2, a shifted-and-inverted spectral transformation replaces the original pencil by a transformed standard eigenproblem. With a shift $A,B\in\mathbb{R}^{d\times d}$3 such that $A,B\in\mathbb{R}^{d\times d}$4 is nonsingular, one factors

$A,B\in\mathbb{R}^{d\times d}$5

forms

$A,B\in\mathbb{R}^{d\times d}$6

and solves

$A,B\in\mathbb{R}^{d\times d}$7

Under assumptions on the shift, the analysis conditionally ensures that every computed shifted and inverted eigenvalue is close to the exact shifted and inverted eigenvalue of a pair of matrices close to $A,B\in\mathbb{R}^{d\times d}$8 and $A,B\in\mathbb{R}^{d\times d}$9. The same analysis yields error bounds for the original problem, and the computed generalized eigenvectors satisfy relative residual bounds that depend on the distance between the eigenvalue and the shift (Stewart, 2024).

Shift choice is therefore structural, not cosmetic. If $$A\phi_i=\lambda_i B\phi_i,\qquad i\in\{1,\dots,d\},$$0 is too close to an eigenvalue, $$A\phi_i=\lambda_i B\phi_i,\qquad i\in\{1,\dots,d\},$$1 becomes nearly singular and conditioning deteriorates. If the shift is of moderate size, the paper states that every computed generalized eigenvalue can correspond to a generalized eigenvalue of a pair of matrices close to $$A\phi_i=\lambda_i B\phi_i,\qquad i\in\{1,\dots,d\},$$2 and $$A\phi_i=\lambda_i B\phi_i,\qquad i\in\{1,\dots,d\},$$3, while relative residuals are small for generalized eigenvalues that are not much larger than the shift (Stewart, 2024).

Singular pencils require a different stabilization logic. A rank-completing perturbation regularizes a singular square pencil $$A\phi_i=\lambda_i B\phi_i,\qquad i\in\{1,\dots,d\},$$4 by adding

$$A\phi_i=\lambda_i B\phi_i,\qquad i\in\{1,\dots,d\},$$5

where $$A\phi_i=\lambda_i B\phi_i,\qquad i\in\{1,\dots,d\},$$6, $$A\phi_i=\lambda_i B\phi_i,\qquad i\in\{1,\dots,d\},$$7 and $$A\phi_i=\lambda_i B\phi_i,\qquad i\in\{1,\dots,d\},$$8 are $$A\phi_i=\lambda_i B\phi_i,\qquad i\in\{1,\dots,d\},$$9 matrices with orthonormal columns, $B=I$0 is regular, and $B=I$1. For true eigenvalues of the original pencil preserved under the perturbation, the corresponding right and left eigenvectors satisfy $B=I$2 and $B=I$3, which allows the true eigenvalues to be separated from the prescribed and random eigenvalues introduced by the perturbation (Hochstenbach et al., 2018). This is presented as a fast and robust alternative to staircase methods.

Spectral certification can also be done without solving the problem exactly. Euclidean-metric Gerschgorin-type inclusion sets for $B=I$4 place every eigenvalue in explicitly computable disks. If the $B=I$5th row of $B=I$6 is strictly diagonally dominant,

$B=I$7

where $B=I$8 and $B=I$9. These disks are easier to compute than earlier chordal-metric sets and support forward error analysis for a computed eigenvalue of a diagonalizable pencil (Nakatsukasa, 2010).

Perturbation theory sharpens the distinction between Hermitian and non-Hermitian cases. For generalized Hermitian eigenvalue problems, quadratic perturbation bounds are proportional to the square of the perturbation norm divided by the spectral gap, while the non-Hermitian case generally does not admit such quadratic stabilization (Nakatsukasa, 2010). One practical implication is that “stable” computed generalized eigenvalues in Hermitian settings can legitimately mean second-order insensitivity to small perturbations, whereas in non-Hermitian settings the same term is necessarily weaker.

4. Iterative and optimization-based eigensolvers

Subspace methods provide another layer of GEVS by enriching the search space with directions that improve robustness. For large sparse symmetric GEPs with $$\underset{\phi}{\text{maximize}}\ \phi^\top A\phi \quad\text{subject to}\quad \phi^\top B\phi=1,$$0 positive definite, the Chebyshev-RQI Subspace (CRS) method combines Chebyshev polynomial filtering with inexact Rayleigh quotient iteration. At iteration $$\underset{\phi}{\text{maximize}}\ \phi^\top A\phi \quad\text{subject to}\quad \phi^\top B\phi=1,$$1, one forms

$$\underset{\phi}{\text{maximize}}\ \phi^\top A\phi \quad\text{subject to}\quad \phi^\top B\phi=1,$$2

and enlarges the subspace according to

$$\underset{\phi}{\text{maximize}}\ \phi^\top A\phi \quad\text{subject to}\quad \phi^\top B\phi=1,$$3

The reduced GEP is then solved by Rayleigh-Ritz extraction (Wang et al., 2022).

The reported numerical behavior is explicitly stabilization-oriented. Across all tested problem sizes, up to $$\underset{\phi}{\text{maximize}}\ \phi^\top A\phi \quad\text{subject to}\quad \phi^\top B\phi=1,$$4, CRS is described as more stable and reliable than Chebyshev-Davidson. For the largest problem, the table in the summary reports that for $$\underset{\phi}{\text{maximize}}\ \phi^\top A\phi \quad\text{subject to}\quad \phi^\top B\phi=1,$$5 eigenpairs Chebyshev-Davidson required $$\underset{\phi}{\text{maximize}}\ \phi^\top A\phi \quad\text{subject to}\quad \phi^\top B\phi=1,$$6 outer iterations, $$\underset{\phi}{\text{maximize}}\ \phi^\top A\phi \quad\text{subject to}\quad \phi^\top B\phi=1,$$7 matrix-vector products, and $$\underset{\phi}{\text{maximize}}\ \phi^\top A\phi \quad\text{subject to}\quad \phi^\top B\phi=1,$$8 seconds, whereas CRS required $$\underset{\phi}{\text{maximize}}\ \phi^\top A\phi \quad\text{subject to}\quad \phi^\top B\phi=1,$$9, $$\underset{\Phi}{\text{maximize}}\ \mathrm{tr}(\Phi^\top A\Phi) \quad\text{subject to}\quad \Phi^\top B\Phi=I.$$0, and $$\underset{\Phi}{\text{maximize}}\ \mathrm{tr}(\Phi^\top A\Phi) \quad\text{subject to}\quad \Phi^\top B\Phi=I.$$1 seconds (Wang et al., 2022). The rationale given is that the Chebyshev filter damps unwanted spectral components while the inexact RQI direction corrects local misalignment.

The optimization literature extends this perspective. In the difference-based formulation, classical first-order methods, the power method, and preconditioned mirror descent are interpreted through surrogate minimization, and an accelerated Split-Merge algorithm is developed for the general GEP setting. Its update uses richer second-order information and takes the form

$$\underset{\Phi}{\text{maximize}}\ \mathrm{tr}(\Phi^\top A\Phi) \quad\text{subject to}\quad \Phi^\top B\Phi=I.$$2

with data-dependent coefficients $$\underset{\Phi}{\text{maximize}}\ \mathrm{tr}(\Phi^\top A\Phi) \quad\text{subject to}\quad \Phi^\top B\Phi=I.$$3 derived from local curvature (Liu et al., 3 Jul 2025). The same work reports up to 4-fold speedup over the power method, and when iterative inexact solvers for $$\underset{\Phi}{\text{maximize}}\ \mathrm{tr}(\Phi^\top A\Phi) \quad\text{subject to}\quad \Phi^\top B\Phi=I.$$4 are used, up to 17x reduction in time compared to Cholesky-based solvers. It also notes that Lanczos can stagnate or fail with clustered eigenvalues, whereas Split-Merge succeeded on all trials in the reported experiments (Liu et al., 3 Jul 2025).

These algorithmic strands share a narrow interpretation of stabilization: not arbitrary damping, but a controlled improvement of convergence behavior, conditioning, and numerical reliability while retaining the target eigenspace.

5. Discretization-driven GEVS in computational mechanics and physics

In immersed explicit dynamics, GEVS is applied directly to element mass matrices of cut elements. The difficulty is that poorly cut elements make the local mass matrix nearly singular, creating spectral outliers in $$\underset{\Phi}{\text{maximize}}\ \mathrm{tr}(\Phi^\top A\Phi) \quad\text{subject to}\quad \Phi^\top B\Phi=I.$$5 and reducing the critical time step

$$\underset{\Phi}{\text{maximize}}\ \mathrm{tr}(\Phi^\top A\Phi) \quad\text{subject to}\quad \Phi^\top B\Phi=I.$$6

The generalized eigenvalue stabilization strategy operates on each cut element $$\underset{\Phi}{\text{maximize}}\ \mathrm{tr}(\Phi^\top A\Phi) \quad\text{subject to}\quad \Phi^\top B\Phi=I.$$7 by solving the local generalized eigenproblem

$$\underset{\Phi}{\text{maximize}}\ \mathrm{tr}(\Phi^\top A\Phi) \quad\text{subject to}\quad \Phi^\top B\Phi=I.$$8

identifying the set $$\underset{\Phi}{\text{maximize}}\ \mathrm{tr}(\Phi^\top A\Phi) \quad\text{subject to}\quad \Phi^\top B\Phi=I.$$9 of bad eigenpairs with $R(A,B;x)=\frac{x^\top Ax}{x^\top Bx}$0, and modifying the mass matrix according to

$R(A,B;x)=\frac{x^\top Ax}{x^\top Bx}$1

The paper states that this shifts each outlier eigenvalue $R(A,B;x)=\frac{x^\top Ax}{x^\top Bx}$2 to $R(A,B;x)=\frac{x^\top Ax}{x^\top Bx}$3 without modifying the rest of the spectrum, and that the resulting stabilization achieves optimal convergence rates and recovers critical time step sizes of equivalent boundary-conforming discretizations, including with weakly enforced Dirichlet boundary conditions by Nitsche’s method or penalty formulations (Bürchner et al., 9 Sep 2025).

The same application area also contains an earlier element-level eigenvalue stabilization (EVS) technique. There, one diagonalizes the cut-element mass matrix $R(A,B;x)=\frac{x^\top Ax}{x^\top Bx}$4, classifies spurious modes via

$R(A,B;x)=\frac{x^\top Ax}{x^\top Bx}$5

and constructs a stabilization matrix

$R(A,B;x)=\frac{x^\top Ax}{x^\top Bx}$6

The summary reports increases of the critical timestep by factors of 2.5 for moderate cuts and by more than 30 for extreme cuts, while maintaining convergence behavior comparable to the reference geometry-fitted spectral element method (Eisenträger et al., 2023). The later GEVS formulation generalizes this selective spectral modification and combines it with the finite cell method to guarantee definiteness of the system matrices (Bürchner et al., 9 Sep 2025).

In domain decomposition for the Biot model, generalized eigenproblems are used to construct coarse spaces that target global low-energy error modes. On each overlapping subdomain $R(A,B;x)=\frac{x^\top Ax}{x^\top Bx}$7, the GenEO method solves

$R(A,B;x)=\frac{x^\top Ax}{x^\top Bx}$8

retains eigenvectors with $R(A,B;x)=\frac{x^\top Ax}{x^\top Bx}$9, and assembles the global coarse space

$A\Phi=B\Phi\Lambda$00

The same paper also uses a pressure stabilization term

$A\Phi=B\Phi\Lambda$01

to ensure stability for piecewise constant pressure discretization. The reported GMRES iteration counts grow moderately as the number of subdomains increases from $A\Phi=B\Phi\Lambda$02 to $A\Phi=B\Phi\Lambda$03, and remain modest across wide ranges of Poisson ratio and permeability, supporting claims of scalability and parametric robustness (Lee, 2022).

Lattice gauge theory provides a different but closely related interpretation of GEVS. The generalized eigenvalue problem for correlation matrices is stabilized not by modifying the pencil itself, but by careful choice of reference time $A\Phi=B\Phi\Lambda$04, by pruning nearly linearly dependent interpolators, by projecting out operator combinations corresponding to small eigenvalues, and by checking the stability of physical outputs under basis reduction. The paper recommends monitoring plateau behavior, varying $A\Phi=B\Phi\Lambda$05, and using bootstrap or jackknife error analysis (0808.1017). In that context, stabilization means reliable extraction of energies and matrix elements from a noisy, nearly singular generalized eigenproblem.

6. Quantum formulations

Quantum algorithms stabilize GEPs by transforming or reparameterizing them so that the physically relevant spectrum can be accessed by quantum subroutines. For Hermitian $A\Phi=B\Phi\Lambda$06 and positive-definite $A\Phi=B\Phi\Lambda$07, quantum phase estimation applies to the transformed Hermitian operator

$A\Phi=B\Phi\Lambda$08

The key efficiency criterion is sparsity of $A\Phi=B\Phi\Lambda$09: if $A\Phi=B\Phi\Lambda$10 is diagonal or block-diagonal, $A\Phi=B\Phi\Lambda$11 preserves the sparse structure, and efficient Hamiltonian simulation becomes possible. The paper identifies Sturm–Liouville problems and linear ideal magnetohydrodynamics as paradigmatic examples, with the latter interpreted as a stability problem for magnetically confined plasmas (Parker et al., 2020).

Variational quantum generalized eigensolvers keep the $A\Phi=B\Phi\Lambda$12-weighted structure explicit rather than eliminating it. The basic loss is the generalized Rayleigh quotient

$A\Phi=B\Phi\Lambda$13

augmented for higher states by $A\Phi=B\Phi\Lambda$14-orthogonality penalties. A full quantum generalized eigensolver then performs quantum gradient descent with

$A\Phi=B\Phi\Lambda$15

The reported 2-qubit simulations successfully recover the generalized eigenvalues of the matrix pencil $A\Phi=B\Phi\Lambda$16, and the numerical experiment indicates that FQGE is robust under Gaussian noise (Liang et al., 2021).

A separate variational approach reduces each single-qubit parameter update to a $A\Phi=B\Phi\Lambda$17 GEP,

$A\Phi=B\Phi\Lambda$18

where the generalized Rayleigh quotient is optimized analytically with respect to a unit-quaternion parameterization of the gate. The paper states that if the denominator matrix becomes non-positive definite because of shot noise, a small positive definite shift can be added to mitigate numerical instabilities. It also reports a bonus finding: a problem having a real-valued solution can be solved more effectively using quantum gates generating a complex-valued state vector (Sato et al., 2023).

In these quantum settings, stabilization is primarily structural. The $A\Phi=B\Phi\Lambda$19-weighted inner product is preserved, the problematic inversion of $A\Phi=B\Phi\Lambda$20 is avoided or controlled, and noise sensitivity is managed through small-dimensional GEP subproblems or sparse Hermitian transformations.

7. Terminological scope and related stabilization phenomena

The term “stabilization” is not used uniformly across all spectral literature. In nonlinear $A\Phi=B\Phi\Lambda$21-Laplacian theory, for example, stabilization refers not to regularization of a matrix pencil but to a threshold effect in existence theory. For boundary value problems of the form

$A\Phi=B\Phi\Lambda$22

with $A\Phi=B\Phi\Lambda$23 continuous, nonnegative, and nondecreasing, there exists a threshold

$A\Phi=B\Phi\Lambda$24

such that for all $A\Phi=B\Phi\Lambda$25 a positive solution exists (Carrillo-Reina et al., 27 Apr 2026). The mechanism is derived from an explicit maximum principle depending on $A\Phi=B\Phi\Lambda$26 and the slab diameter $A\Phi=B\Phi\Lambda$27, and from the fixed-point map

$A\Phi=B\Phi\Lambda$28

The same paper gives a new lower bound for the first Dirichlet eigenvalue,

$A\Phi=B\Phi\Lambda$29

and interprets the large-$A\Phi=B\Phi\Lambda$30 behavior as suggesting a connection to the $A\Phi=B\Phi\Lambda$31-Laplacian (Carrillo-Reina et al., 27 Apr 2026).

This suggests that “stabilization” has become a broader term for threshold phenomena in spectral or eigenvalue-controlled problems. In matrix-based GEVS, however, the core meaning remains narrower: stabilization is the deliberate control of generalized eigenvalues, eigenvectors, or the weighting operator so that the computation remains meaningful under ill-conditioning, singularity, discretization artifacts, or noise. Under that narrower meaning, the literature consistently emphasizes selective intervention—regularizing $A\Phi=B\Phi\Lambda$32, choosing a shift, completing rank, enriching a subspace, or modifying only the bad local modes—rather than indiscriminate spectral damping.