Polynomial Eigenvalue Method

• Polynomial eigenvalue method is an algorithmic framework that finds eigenvalues and eigenvectors of matrix polynomials without explicit determinant computation.
• It employs trace-based Laguerre iterations and QR factorizations to achieve reliable, backward-stable convergence in applications like vibration analysis and control theory.
• The method exploits structural sparsity (e.g., Hessenberg and tridiagonal cases) to reduce computational complexity while maintaining high numerical accuracy.

A polynomial eigenvalue method refers to algorithmic frameworks for determining the eigenvalues and corresponding eigenvectors of matrix polynomials—expressions of the form $P(\lambda) = \sum_{i=0}^d \lambda^i A_i$ with $A_i \in \mathbb{C}^{n \times n}$ and $A_d \neq 0$. The central problem involves finding scalars $\lambda$ such that $\det P(\lambda) = 0$ and nonzero vectors $x$ satisfying $P(\lambda)x = 0$. This problem arises in vibration analysis, control theory, and diverse engineering applications. The eigenvalues may be finite or infinite, with the latter associated to the reversal polynomial. Modern polynomial eigenvalue methods eschew forming the determinant explicitly, instead relying on stable iterative techniques—most notably the Laguerre root-finding iteration adapted for matrix polynomials (Cameron et al., 2017).

1. Formulation and Structure of the Matrix Polynomial Eigenvalue Problem

The matrix polynomial eigenvalue problem (PEP) is formalized as seeking roots $\lambda$ of the scalar polynomial $p(\lambda)=\det P(\lambda)$, which encompasses both finite and infinite eigenvalues. Infinite eigenvalues correspond to zeros of the reversal polynomial, $$\mathrm{rev}\,P(\rho) = \sum_{i=0}^d \rho^{d-i} A_i$$, with $\rho=1/\lambda$. The eigensystem encompasses $nd$ eigenvalues (including multiplicities), and associated eigenvectors $x$ such that $P(\lambda)x=0$. Particularly, the underlying algebraic multiplicity is governed by the regularity condition $A_d\neq 0$.

Traditional approaches proceed via linearization, embedding $P(\lambda)$ in a larger matrix pencil $L(\lambda) = A - \lambda B$ to exploit well-established GEP solvers. However, such linearization can compromise numerical conditioning and stability in practice (You et al., 2017). Purpose-built iterative root-finding methods act directly on the matrix polynomial's spectral structure, enabling backward-stable computation with quantifiable error control (Cameron et al., 2017, Bini et al., 2012).

2. Laguerre Iteration Applied to Matrix Polynomials

Distinctive among root-finding techniques for matrix polynomials is the trace-based Laguerre iteration. The scalar Laguerre method employs

$$S_1(\lambda) = \frac{p'(\lambda)}{p(\lambda)}, \quad S_2(\lambda) = -\bigl(S_1(\lambda)\bigr)' = \sum_{i=1}^{N_1} \frac{1}{(\lambda - r_i)^2}$$

where the update step is

$$\hat{\lambda} = \lambda - \frac{N_1}{S_1(\lambda) \pm \sqrt{ (N_1-1)(N_1 S_2(\lambda) - S_1(\lambda)^2 )}}$$

In the matrix setting, $S_1(\lambda)$ and $S_2(\lambda)$ are stably computed via Jacobi's formula and matrix equations:

$$P(\lambda) X_1(\lambda) = P'(\lambda), \quad P(\lambda) X_2(\lambda) = P''(\lambda)$$

$$S_1(\lambda) = \operatorname{trace}(X_1(\lambda)),\quad S_2(\lambda) = \operatorname{trace}( X_1(\lambda)^2 - X_2(\lambda) )$$

Each iteration thus solves two $n\times n$ linear matrix equations—most conveniently via QR factorization—and evaluates only their trace, avoiding direct computation of $\det P(\lambda)$. The iteration proceeds until convergence is detected by residual norm, backward-error bound, or stagnation in the update (Cameron et al., 2017).

3. Initialization via Numerical Range and Newton Polygon

Robust convergence of any simultaneous root-finding method critically depends on high-quality initial guesses for eigenvalues. Two principal approaches are detailed:

• Newton-polygon radii: Construct the scalar polynomial $w(\lambda)=\sum_{i=0}^d \|A_i\|\,\lambda^i$. The Newton polygon yields a sequence of radii; abscissas $0 = k_0 < k_1 < \dots < k_q = d$ supply points on circles $|\lambda|=r_i$ where $$r_i = |\frac{a_{k_{i-1}}}{a_{k_i}}|^{1/(k_i-k_{i-1})}$$, for $a_j = \|A_j\|$.
• Numerical-range roots: Define $W(P) = \{\lambda \in \mathbb{C} : x^*P(\lambda)x=0$ for some $\|x\|=1\}$. Using column-pivoted QR of $A_0$ and $A_d$, scalar polynomials in the form $q_j^*P(\lambda)q_j = 0$ are solved for each column $q_j$, producing root clusters aligned to spectral features of $P$. For self-adjoint $P$, eigenvalues in $W(P)$ are bounded by upper Pellet bounds.

This principled selection of starting points is especially critical for highly clustered spectra and underpins the rapid global convergence witnessed in practice (Cameron et al., 2017).

4. Stopping Criteria, Convergence, and Backward Stability

Iteration on candidate eigenvalues is halted when one of three criteria is met:

• Small residual: If the smallest diagonal of the QR factorization of $P(\lambda)$ falls below $\tau$, convergence is declared. The threshold is set by

$$\tau = \begin{cases} \alpha\,\epsilon, & |\lambda| \leq 1 \ \tilde{\alpha}\,\epsilon, & |\lambda| > 1 \end{cases}$$

with $\alpha = \sum_{i=0}^d |\lambda|^i \|A_i\|$, $\tilde{\alpha} = \sum_{i=0}^d |\rho|^{d-i} \|A_i\|$.

• Backward-error bound: If $$\min_{b\neq 0} \|b\|\, \alpha \|P(\lambda)^{-1} b\| < \epsilon$$ (the Tisseur bound), $\lambda$ is accurate to precision $\epsilon$.
• Stagnation: If $|\hat{\lambda} - \lambda| < \epsilon |\lambda|$, or iteration reaches its maximum.

For computed $\lambda$ meeting either of the first two, the eigenpair is guaranteed backward stable to within $O(n\epsilon)$. This explicitly quantifies the numerical reliability of the computed spectrum (Cameron et al., 2017).

5. Computation of Eigenvectors and Condition Numbers

Upon convergence of $\lambda$:

• Right and left eigenvectors are constructed using the QR factorization. If a small pivot exists, the right eigenvector $x$ is obtained by solving a reduced system, with the left eigenvector $y$ set to the final column of $Q$. Otherwise, inverse iteration is applied to $R^*R$ or $RR^*$.
• The normwise condition number of the simple eigenpair $(\lambda,x,y)$ is

$$\kappa(\lambda, P) = \frac{\alpha \|x\|\|y\|}{|\lambda y^*P'(\lambda)x|}$$

with analogous expressions for zero or infinite eigenvalues. This facilitates sensitivity analysis and quantifies perturbation effects (Cameron et al., 2017).

6. Exploiting Structural Sparsity: Hessenberg and Tridiagonal Cases

Significant computational savings arise when the sequence $\{A_i\}$ is structured (Hessenberg or tridiagonal):

• Hyman's method evaluates $\det P(\lambda)$ and derivatives in $O(n^2)$ (Hessenberg) or $O(n)$ (tridiagonal), rather than $O(n^3)$ per evaluation.
• The approach analyzes the last row of $P(\lambda)$, employing forward and adjoint solves to assemble determinant and derivative quantities, delivering log-derivative formulas:

$$\frac{p'(\lambda)}{p(\lambda)} = \frac{b'(\lambda)}{b(\lambda)} + \frac{q'(\lambda)}{q(\lambda)}$$

This reduces per-iteration cost sharply relative to general dense cases, with tridiagonal complexity scaling as $O(d n^2 + d^2 n^2)$. Efficient core operations enable the method’s viability for large-scale or highly structured eigenproblems typical in applied settings (Cameron et al., 2017).

7. Performance, Numerical Experiments, and Scope

Benchmarks confirm:

• For scalar polynomials ($n=1$), the method matches or outperforms established algorithms (POLZEROS, AMVW) in speed and accuracy, with $O(d^2)$ scaling.
• For tridiagonal and Hessenberg problems, computational cost corresponds to $O(d^2 n^2)$ and $O(dn^3 + d^2 n^3)$, respectively—substantially reduced from dense cases.
• Backward errors remain within a few units of machine precision, matching premier solvers (QUADEIG).
• Detailed forward error plots and tables demonstrate stability and competitiveness across a wide spectrum of eigenvalue problem sizes and structures.

The combination of Laguerre’s method, trace-based evaluation, spectrum-respecting initialization, and efficient structure exploitation yields a unified, backward-stable, and robust computational apparatus for polynomial eigenvalue problems (Cameron et al., 2017).

---

References

• Cameron, Steckley. "On the application of Laguerre's method to the polynomial eigenvalue problem" (Cameron et al., 2017).
• Bini, Noferini. "Numerical methods for polynomial root-finding."
• Tisseur. "Backward error analysis in matrix polynomials."
• Plestenjak, Bini et al., "Tridiagonal and Hessenberg solvers for matrix polynomials."
• AMVW: Aurentz, Mach, Vandebril, Watkins. "Structured Polynomial Eigenvalue Problems."
• QUADEIG: Aurentz, Mach, Vandebril, Watkins, "Fast and backward stable computation of the eigenvalues and eigenvectors of matrix polynomials."