Nonlinear Matrix Eigenvalue Problem

Updated 13 December 2025

Nonlinear matrix eigenvalue problems involve finding nontrivial solutions (x, λ) such that T(λ)x = 0, where T(λ) is a holomorphic matrix-valued function with nonlinear dependence on λ.
They arise in diverse applications, including quantum mechanics, photonics, delay-differential equations, and control theory, where advanced spectral analysis is critical.
Modern numerical approaches utilize contour-integral methods, Krylov subspace techniques, and rational approximations to accurately localize eigenvalues and ensure robust convergence.

A nonlinear matrix eigenvalue problem (NMEP) is the task of computing nontrivial vectors $x \in \mathbb{C}^n$ and associated scalars $\lambda \in \mathbb{C}$ such that $T(\lambda)x = 0$ , where $T: \Omega \to \mathbb{C}^{n\times n}$ is a holomorphic (often analytic, rational, or meromorphic) matrix-valued function defined on some domain $\Omega \subset \mathbb{C}$ . NMEPs generalize the linear matrix eigenproblem by permitting arbitrary (not necessarily linear or polynomial) dependence on the spectral parameter $\lambda$ , and arise in a wide range of applications including quantum mechanics, photonics, delay-differential equations, structural vibrations, and control theory. Their mathematical and computational analysis combines complex operator theory, nonlinear approximation, numerical linear algebra, and high-performance algorithms.

1. Mathematical Foundations

The canonical NMEP is

$T(\lambda)x = 0, \quad x \neq 0,$

where typically $T(\lambda)$ is analytic on an open subset $\Omega \subset \mathbb{C}$ . The spectrum $\Lambda(T) = \{ \lambda \in \Omega \mid \det T(\lambda) = 0 \}$ consists of isolated points (possibly accumulating at the boundary or at infinity), each counted with algebraic multiplicity determined by the order of vanishing of $\det T$ at $\lambda$ . Keldysh's theorem extends the structure theory of the spectrum: the matrix function $T(\lambda)^{-1}$ , viewed as a meromorphic function in the resolvent set, admits a partial fraction expansion whose singular terms encode the right and left eigenvectors and Jordan structure of $T$ (Beyn, 2010).

In the context of infinite-dimensional spaces (e.g., operator-valued NMEPs representing PDEs or delay equations), $T(\lambda)$ becomes a family of closed, Fredholm operators acting on a Hilbert space, and the discrete eigenvalues correspond to isolated singularities of $T(\lambda)^{-1}$ . The pseudospectrum, and related localization theorems, generalize to the nonlinear setting, and play a crucial role in spectral verification and error analysis (Bindel et al., 2013, Colbrook et al., 2023).

2. Classification and Model Problems

NMEPs naturally stratify according to the analytic structure of $T(\lambda)$ and the source of nonlinearity:

Polynomial eigenvalue problems (PEPs): $T(\lambda) = \sum_{k=0}^d \lambda^k A_k$ , which admit companion linearizations and are a classical subject.
Rational eigenvalue problems (REPs): $T(\lambda) = P(\lambda) + \sum_{i=1}^m (C_i - \lambda D_i) r_i(\lambda)$ , with $r_i$ rational (Lietaert et al., 2018, Saad et al., 2019, Aziz et al., 2020).
Transcendental-type NMEPs: $T(\lambda)$ contains terms such as exponentials or roots, e.g., $e^{-\sqrt{\lambda}}$ (Lietaert et al., 2018).
Eigenvector-nonlinear NMEPs (NEPv): The matrix depends explicitly on the eigenvector, $T(x, \lambda)x = 0$ (Jarlebring et al., 2012, Jarlebring et al., 2020, Janssens et al., 3 Oct 2025).
Infinite-dimensional/Operator NMEPs: $T(\lambda): D(T)\rightarrow H$ , $H$ Hilbert, often arising from spatially discretized or semi-discretized PDEs (Colbrook et al., 2023).

Many physically important cases reduce to block or Schur complements of more general NMEPs, and some nonlinearities enter via domain decomposition or parametric model-reduction schemes (Ringh et al., 2019).

3. Analytical and Perturbation Theory

The spectral and perturbation theory for NMEPs is more intricate than in the linear case. The nonlinear generalization of the Bauer–Fike theorem provides a sharp bound for eigenvalue perturbations: $| \tilde{\lambda} - \lambda_0 | \leq \kappa(T, \lambda_0) \cdot \| T(\lambda_0) - \tilde{T}(\tilde{\lambda}) \| + O(\|T-\tilde{T}\|^2),$ where $\kappa(T, \lambda_0) = \|T'(\lambda_0)^{-1}\|$ is a nonlinear condition number (Katende, 17 Sep 2024, Bindel et al., 2013). Keldysh’s theorem justifies representing $T(z)^{-1}$ as a Laurent expansion near simple or semi-simple eigenvalues, enabling the definition of spectral projectors and invariance subspaces via contour integration.

For eigenvector-nonlinear NMEPs, the appropriate notion of conditioning and convergence is subtler. Local convergence rates for iterative methods can often be derived from the linearization of the map $x\mapsto T(x, \lambda)x$ at a solution, and closely mirror those of the associated Jacobian (Jarlebring et al., 2012, Jarlebring et al., 2020).

Localization theorems generalizing Gershgorin’s disks and the pseudospectral inclusion for analytically parameterized matrix functions enable a priori demarcation of regions containing eigenvalues and rigorous eigenvalue counting (Bindel et al., 2013).

4. Numerical Methods

A variety of algorithms have been developed for NMEPs, each targeting different structures and scales. The primary classes are:

Contour-integral and moment-based methods:

Beyn's integral method computes all eigenvalues inside a contour $\Gamma$ by contour integrating resolvent-moment matrices, reducing the problem to a small projected linear eigenproblem whose spectrum matches the targeted eigenvalues (Beyn, 2010, Xi et al., 2023, Brennan et al., 2020).
Riesz-projection-based methods use scalar functionals of the resolvent to form moment equations for the eigenvalues within a contour, allowing prioritization of physically relevant (coupled) modes without linearization (Binkowski et al., 2018).

Krylov–Arnoldi–Lanczos extensions:

Infinite bi-Lanczos and variants construct Krylov subspaces for infinite-dimensional linearizations of $T(\lambda)$ to simultaneously approximate right and left eigenvectors (Gaaf et al., 2016).
Rational/truncated Krylov methods (e.g., CORK, TS-CORK) exploit rational approximations and compact linearizations, achieving efficient shift-and-invert strategies for large-scale REPs (Lietaert et al., 2018, Aziz et al., 2020).

Quasi-Newton, Newton, and Broyden-type methods:

Newton-type iteration for $(\lambda, x)$ seeks stationary points of the system $F(\lambda, x) = (T(\lambda)x, c^H x - 1)$ , using Block-Newton or quasi-Newton linearizations in each iteration (Jarlebring et al., 2017). Several classical methods can be interpreted as quasi-Newton schemes, such as residual inverse iteration and the method of successive linear problems, each with quantifiable convergence rates based on Keldysh-theoretic resolvent expansions.
Broyden's method adapts the secant-based approach to the structured augmented NMEP system, achieving local superlinear convergence for simple eigenvalues and supporting robust deflation strategies (Jarlebring, 2018).

Rational approximation and linearization:

AAA rational approximations and set-valued adaptations construct rational surrogates matching the nonlinear terms in $T(\lambda)$ , enabling efficient CORK-style companion linearizations and resulting in pencils of minimal dimension that embed the true spectrum (Lietaert et al., 2018, Saad et al., 2019).

Multigrid and multilevel schemes:

Multilevel correction methods decompose the nonlinear eigenproblem into a sequence of linear boundary-value problems on multigrid hierarchies and small-dimensional nonlinear corrections, reducing overall computational work to that of a linear problem of the same scale (Jia et al., 2015, Xie, 2014).

NEPv-focused algorithms:

For eigenvector-dependent problems, implicit Newton schemes and structure-exploiting inverse iterations form the prevailing methodology. Quadratic convergence is achievable for Newton-type updates (proper Jacobian handling) while self-consistent field (SCF) iterations offer linear, robust progress (Jarlebring et al., 2012, Jarlebring et al., 2020).

A summary table of selected method classes, key ideas, and cost structure:

Methodology	Core Principle	Complexity per Iteration
Contour methods	Resolvent moment integrals	$O(N \cdot l)$ solves, $N \sim 30$ –$50$
Rational Krylov	Linearization via AAA/CORK	$O(n (\text{block size})^2)$ per step
Quasi-Newton	Fixed-shift Jacobian solve	1 LU + 1 mat-vec or eigenproblem
Multigrid FMG	Nested mesh linear solves	$O(\mathrm{total \; dof})$ per full solve
Eigenvector NEP	Eigendecomposition of Jacobian at $v$	$O(n^3)$

5. Implementation Issues and Spectral Verification

Critical challenges arise in avoiding artifacts due to discretization of operator-valued NMEPs:

Spectral pollution: Discretized spectra may exhibit spurious modes, poor convergence, or missing eigenvalues, especially when dealing with unbounded operators or non-compact resolvents (Colbrook et al., 2023). The infBeyn method circumvents this by constructing the contour-integral reduced pencil directly in the operator setting and only discretizing after projection.
Pseudospectra and stability: Localizing eigenvalues robustly is achievable via nonlinear analogues of the pseudospectrum, featured in both localization theorems and the explicit error control in infinite-dimensional settings (Bindel et al., 2013, Colbrook et al., 2023).
Deflation and multiple eigenvalues: Augmented or block-deflated systems permit systematic computation of several (or all) modes, given suitable orthogonality conditions (Jarlebring, 2018, Beyn, 2010).

6. Practical Applications and Example Problems

NMEPs are central in several domains:

Quantum mechanics: Nonselfadjoint boundary value problems (e.g., resonances in open quantum wells or nanophotonic structures) frequently reduce to rational or transcendental NMEPs (Binkowski et al., 2018).
Photonic and nanophotonic systems: Maxwell or Helmholtz equations with dispersive, frequency-dependent permittivities yield large-scale rational NMEPs after finite element discretization (Lietaert et al., 2018, Binkowski et al., 2018).
Delay-differential and control systems: Time-delay in feedback leads to matrix-valued functions involving exponentials or other transcendental functions (Gaaf et al., 2016, Xi et al., 2023).
Gross–Pitaevskii and nonlinear Schrödinger equations: NEPv with quadratic or rational dependence, especially in the presence of normalization and symmetry properties, appear in Bose–Einstein condensate modeling (Jarlebring et al., 2012, Jarlebring et al., 2020, Janssens et al., 3 Oct 2025).

Designed numerical experiments confirm that the most advanced methods recover all eigenvalues inside a given region to machine precision, up to sizes $n \sim 10^5$ – $10^6$ , and efficiently separate physical from spurious modes in physics-informed settings (Beyn, 2010, Binkowski et al., 2018, Lietaert et al., 2018, Colbrook et al., 2023).

7. Advanced Topics and Recent Developments

Several frontiers are being actively pursued:

Subspace and derivative-interpolatory frameworks combine model reduction concepts with Hermite interpolation to achieve quadratic (or higher) convergence for select eigenvalues, especially in large-scale rational NMEPs (Aziz et al., 2020).
Multiparameter linearizations for NEPv with rational eigenvector dependencies realize spectral inclusion via structured operator determinants and exploit tailored Krylov–Arnoldi filtering to mitigate the swamping of spurious solutions inherent in generic linearizations (Janssens et al., 3 Oct 2025).
Adaptive and parallel contour-integration: Domain decomposition and parallelized quadrature permit scaling to complex eigenvalue landscapes without a priori enumeration of eigenvalue count or distribution (Xi et al., 2023).
Sensitivity, conditioning, and stability: Recent advances extend the Bauer–Fike theorem explicitly to nonlinear and operator-valued settings, providing rigorous frameworks for step-size control in iterative solvers, verification of computed eigenvalues, and parametric tracking in bifurcation scenarios (Katende, 17 Sep 2024).
Automatic rational approximation: Black-box frameworks (e.g., AAA-CORK) automate pole/zero assignment for linearization while achieving optimal linearization sizes and negligible user overhead (Lietaert et al., 2018).

Open research challenges include effective iterative solvers for massive-scale operator-NMEPs, universal spectral localization tools for highly defective spectra, structure-preserving generalizations for non-Hermitian or indefinite problems, and integration with parameter continuation and bifurcation tracking mechanisms.

References (arXiv IDs):