Nonlinear Eigenvalue Problems: Theory & Methods

• Nonlinear eigenvalue problems are defined by holomorphic matrix functions whose singularities determine eigenvalues, generalizing classical linear problems.
• They are applied in areas such as vibration analysis, photonics, quantum resonance, and electronic structure computations, demonstrating wide practical relevance.
• Advanced numerical strategies like rational linearization, contour integration, and homotopy continuation enable accurate and scalable solutions.

A nonlinear eigenvalue problem (NEP) is defined by a holomorphic matrix-valued function $$T : \Omega \subset \mathbb{C} \rightarrow \mathbb{C}^{n \times n}$$. The NEP seeks values $\lambda \in \Omega$ and nonzero vectors $x \in \mathbb{C}^n$ such that $T(\lambda)x = 0$. Nonlinear eigenvalue problems generalize the linear eigenproblem by admitting general holomorphic, rational, algebraic, or transcendental parameter dependence in $T(\lambda)$. These problems arise naturally in numerous application domains: vibration analysis (quadratic and polynomial EVP), photonics (rational and delay NEP), quantum resonance computation, electronic structure (Kohn–Sham, Hartree–Fock), nonlocal pattern formation, and topological systems with eigenvalue-dependent Hamiltonians. NEPs are structurally richer and more computationally challenging than their linear counterparts, and their solution, localization, and spectral theory have become central to modern computational mathematics and physics.

1. Mathematical Foundations and Spectral Theory

The generic NEP takes the form $T(\lambda)x = 0$, with holomorphic (or at minimum continuous) $T$. Often, $T$ is given in split form: $T(\lambda) = \sum_{i=0}^d A_i f_i(\lambda),$ where the $A_i$ are constant matrices and $\lambda \in \Omega$0 are scalar functions (polynomial, rational, algebraic, etc) (Campos et al., 2019).

The spectrum of $\lambda \in \Omega$1 consists of those $\lambda \in \Omega$2 for which $\lambda \in \Omega$3 is singular, i.e., $\lambda \in \Omega$4, and left eigenvectors $\lambda \in \Omega$5 such that $\lambda \in \Omega$6 may also be of interest. Important theoretical tools include:

• Keldysh’s theorem: For holomorphic $\lambda \in \Omega$7, this establishes contour-integral representations for the spectral projector onto the eigenspaces within a domain, as well as local meromorphic resolvent expansions (Bindel et al., 2013, Colbrook et al., 2023).
• Generalized Gershgorin, pseudospectral, and Bauer–Fike theorems: The spectrum can be inclusively localized by nonlinear generalizations of these classical results. For instance, the spectrum is contained within unions of (possibly disconnected) nonlinear Gershgorin regions parameterized by row and column sums of the analytic off-diagonal perturbation, and pseudospectra provide robust error inclusions under perturbation (Bindel et al., 2013).
• Infinite-dimensional setting: In operator-theoretic NEPs, operators $\lambda \in \Omega$8 may be unbounded but densely defined; the discrete spectrum is then governed by analytic Fredholm theory and appropriate decay/compactness conditions (Colbrook et al., 2023).

Nonlinear eigenvalue multiplicity and the structure of Jordan chains critically influence both theoretical analysis and algorithmic convergence (Ringh et al., 2019).

2. Algorithmic Methodologies

A diverse set of numerical strategies exists for NEPs, tailored by the analytic structure of $\lambda \in \Omega$9, problem size, and the desired spectral target. The principal approaches are as follows:

2.1 Rational and Polynomial Linearization

Polynomial eigenvalue problems (PEP), such as

$x \in \mathbb{C}^n$0

are commonly embedded into a companion-type block linearization, reducing the NEP to a (large) generalized linear eigenvalue problem (GEP), for which established linear solvers (QZ, Krylov–Schur) apply (You et al., 2017, Campos et al., 2019). Rational NEPs admit analogous expansions with auxiliary variables for each pole (Lietaert et al., 2018, Saad et al., 2019).

Homotopy continuation is a direct nonlinear alternative, constructed as a path-tracking method between a start-system with known solutions and the target PEP, certified via Smale’s $x \in \mathbb{C}^n$1-theory for correctness. This yields all eigenpairs with near-optimal backward error and can be more stable and better conditioned than linearization, at the cost of higher (but parallelizable) computational effort (You et al., 2017).

2.2 Contour Integral-Based Solvers

Contour integral algorithms exploit the spectral projector representation: $x \in \mathbb{C}^n$2 to extract all eigenvalues within a prescribed domain $x \in \mathbb{C}^n$3. Principal methods include:

• Beyn’s method: Computes moment matrices from randomized probes and contour integration, then solves a reduced GEP; lacks iterative refinement, making high-accuracy computation quadrature-cost limited (Kressner et al., 11 Jun 2026).
• FEAST/NLFEAST: Iterative approaches embedding a nonlinear Rayleigh–Ritz extraction within a contour-filtering framework. NLFEAST achieves linear convergence: for fixed quadrature nodes $x \in \mathbb{C}^n$4 and contraction rate $x \in \mathbb{C}^n$5, subspace error is reduced by $x \in \mathbb{C}^n$6 per iteration. NLFEAST is “embarrassingly parallel” in both quadrature nodes and spectral intervals, efficiently isolates all eigenpairs in a region, and outperforms non-iterative methods such as Beyn at high accuracy (Gavin et al., 2018, Kressner et al., 11 Jun 2026).
• Riesz-projector approaches: Circumvents large linearizations: eigenvalue information is recovered from scalar contour integrals using weighted observables; only a small dimension-nonlinear equation system needs to be solved to recover physically relevant eigenvalues. This is effective when only a subset of eigenpairs matter (e.g., in resonance filtering) (Binkowski et al., 2018).

2.3 Rational Approximation and Compact Rational Krylov

General nonlinearities can be approximated by adaptive barycentric rational interpolants (AAA algorithm), resulting in a rational NEP which admits an explicit block-linearization (CORK, TS-CORK). The set-valued AAA yields compact, automatically-adapted surrogates, greatly reducing the linearization size for a given accuracy relative to Newton or Leja–Bagby bases (NLEIGS). CORK-based rational Krylov methods exploit this structure for scalable per-iteration cost and are highly efficient for large-scale NEPs (Lietaert et al., 2018).

2.4 Projection and Subspace Methods

Generalizations of Arnoldi, infinite Arnoldi (IAR), and bi-Lanczos methods operate on infinite-dimensional operator representations, with finite representations for practical computation. These methods give simultaneous right and left approximations and are effective for interior eigenpairs (Gaaf et al., 2016). Nonlinear Arnoldi and block variants are implemented in parallel NEP solvers (Campos et al., 2019).

2.5 Newton-Type and Quasi-Newton Methods

Newton iterations for NEPs are formulated on the augmented system $x \in \mathbb{C}^n$7, with a block Jacobian. A hierarchy of quasi-Newton methods is obtained by freezing blocks (QN1, QN2, residual inverse iteration, and method of successive linear problems), trading off per-iteration cost versus local convergence rate. QN2 and residual inverse iteration achieve linear convergence, while the successive linear problem method is quadratic but more expensive per step (Jarlebring et al., 2017, Jarlebring, 2018). Structured Broyden variants further reduce work and allow efficient deflation of computed eigenpairs (Jarlebring, 2018).

2.6 Multigrid and Variational Methods for Nonlinear PDE EVPs

Nonlinear eigenproblems arising from variational formulations (e.g., Kohn–Sham, Hartree–Fock, nonlinear optics, pattern formation) are efficiently solved by multigrid approaches. The “multilevel correction” and “Newton-multigrid” algorithms perform a single nonlinear solve on the coarsest mesh; corrections at finer grids solve linearized boundary-value problems, achieving optimal complexity $x \in \mathbb{C}^n$8 where $x \in \mathbb{C}^n$9 is the dimension of the finest discretization (Xie, 2014, Jia et al., 2015, Xu et al., 2024). Damped Newton–mixing and backtracking ensure global convergence even for strong nonlinearities (Xu et al., 2024, Jia et al., 2015).

3. Localization, Conditioning, and Backward Stability

Comprehensive spectral inclusion regions can be computed using nonlinear generalizations of Gershgorin disks and pseudospectral sets; these theorems support both a priori eigenvalue estimates and a posteriori certification of computed values (Bindel et al., 2013). For PEPs and NEPs, conditioning and backward error warrant careful analysis:

• Direct solution of nonlinear problems (e.g., via homotopy continuation for PEPs) can yield better conditioned and more accurate eigenpairs, especially for heavily damped or nonmonic problems, than linearization-based approaches (You et al., 2017).
• When transforming two-parameter linear EVPs to NEPs (e.g., via domain decomposition), the overall conditioning may be benign, provided the small GEP is well-conditioned and solved accurately (Ringh et al., 2019).
• Infinite-dimensional NEPs require special care: naive discretization (e.g., finite element) can introduce spurious eigenvalues (“spectral pollution”), entirely miss spectral components (“invisibility”), or lead to catastrophic ill-conditioning. Discretization-aware contour algorithms (InfBeyn) and direct computation of pseudospectra with explicit error control provide robust alternatives (Colbrook et al., 2023).

A scaled backward error

$T(\lambda)x = 0$0

generalizes the linear residual measure and is widely used for NEP solver tolerance (Campos et al., 2019).

4. Application Domains and Special Problem Classes

NEPs arise in a wide variety of applications, each suggesting distinct analytic or computational approaches.

• Polynomial (QEP/PEP): Vibrational analysis, structural mechanics, delay differential equations. These are frequently amenable to linearization, homotopy, and rational Krylov methodologies (You et al., 2017).
• Rational/Algebraic: Scattering resonances (Maxwell, Helmholtz), nano-photonics, electromagnetic cavities. Rational approximations, NLEIGS, AAA-CORK, and Riesz projections are central (Saad et al., 2019, Binkowski et al., 2018, Lietaert et al., 2018).
• Delay/Transcendental: Control, viscoelasticity, and milling chatter—often require specialized projection/interpolation or operator-based approaches (Bindel et al., 2013, Jarlebring, 2018).
• Nonlocal nonlinearities: Variational formulations in convolution/operator equations, found in nonlocal pattern formation and some quantum models. Multigrid and iterative improvement are preferred (Herrmann et al., 2019).
• Topological and Boundary-sensitive Systems: In systems with eigenvalue-dependent Hamiltonians, such as nonlinear Chern insulators, traditional Bloch theory fails, and the spectrum is extraordinarily sensitive to boundary conditions. Non-Bloch band theory, through generalized Brillouin zones and nonlinear Chern invariants, restores a bulk–boundary correspondence in the nonlinear setting (Otsuka et al., 27 Apr 2026).

5. Parallel and High-Performance Software Ecosystems

NEP solvers have been implemented in highly scalable, parallel libraries such as SLEPc. Key features:

• Flexible problem definition via split form or user-supplied analytic callbacks for $T(\lambda)x = 0$1.
• Parallel data structures (PETSc Mat/Vec objects).
• Support for Newton-type, residual inverse iteration, nonlinear Arnoldi, polynomial/rational interpolation (Chebyshev, NLEIGS), block Krylov–Schur, and TOAR-based compact rational Krylov methods.
• Automatic deflation, region-based spectral filtering, and robust backward-error reporting.
• Demonstrated scalability up to millions of unknowns and hundreds of MPI ranks, with rational interpolation (NLEIGS/AAA–CORK) consistently outperforming other algorithms in strong scaling tests for large-scale NEPs with rational or transcendental parameter dependence (Campos et al., 2019).

Limitations include handling of complex spectra in real arithmetic, deflation accuracy in block methods, and user intervention required for pole/interpolation degree selection in certain algorithms—though adaptive AAA and Leja–Bagby strategies ameliorate some of these concerns.

6. Recent Developments and Advanced Phenomena

Recent research emphasizes:

• Adaptive, automatic rational approximation (AAA, Set-valued AAA): These strategies minimize the linearization dimension for a given error on mesh-independent domains, outperforming fixed-pole Newton bases (Lietaert et al., 2018).
• Robust infinite-dimensional solvers: InfBeyn, combined with precise pseudospectrum computation, has addressed longstanding issues of spectral pollution and invisibility in operator NEPs (Colbrook et al., 2023).
• Non-Bloch frameworks: The emergence of non-Bloch band theory has revealed novel spectral sensitivity and restored topology in NEPs with eigenvalue-dependent Hamiltonians, extending topological band theory beyond conventional Bloch paradigms (Otsuka et al., 27 Apr 2026).
• Rigorous guarantees: Certification via Smale’s $T(\lambda)x = 0$2-theory and explicit error bounds for both direct homotopy and contour-based solvers has strengthened the reliability of NEP computations, especially where spectral clustering or accumulation is present (You et al., 2017, Colbrook et al., 2023).

7. Outlook and Open Problems

NEPs remain an active and rapidly developing area of numerical analysis and applied mathematics, driven both by novel mathematical challenges and emerging application domains. Outstanding research directions include:

• General theory for the conditioning of highly nonlinear and infinite-dimensional problems, especially in the presence of essential spectrum or branch-point singularities.
• Robust, fully automatic pole and basis selection for rational approximation methods in arbitrary domains.
• Unified frameworks for nonlocal, nonlinear, and non-Bloch operator problems with guaranteed spectral inclusivity and a posteriori certification.
• Extension and integration of operator-theoretic perspectives (e.g., Keldysh theory, infinite operator pencils) with efficient high-performance computing infrastructures.

Together, these developments have expanded both the breadth of solvable NEPs and the rigor with which their spectra can be accurately located, computed, and interpreted across the physical and mathematical sciences.