Source-linked AI summary

An integral method for solving nonlinear eigenvalue problems

Wolf-Jürgen Beyn

arXiv:1003.1580v1math.NA

TL;DR

The paper tackles global computation of eigenvalues and eigenvectors for nonlinear holomorphic eigenproblems inside a contour. It uses resolvent contour integrals and Keldysh’s theorem to reduce the problem to a smaller linear eigenproblem, with extensions for difficult multiplicity cases. Numerical examples show that about 25 quadrature nodes often provide good estimates, while the method has a moderately-large-problem scope and can fail in rank-defective cases.

  • Problem

    Existing approaches require eigenvalue localization and iterative computation, while polynomial truncation depends on convergence properties and Newton-type methods remain local.

  • Method

    The method evaluates resolvent contour integrals applied to a random k-column matrix and uses Keldysh’s theorem to reduce the problem to a linear eigenproblem of dimension k.

  • Results

    N ≈25 quadrature nodes usually suffice for good eigenvalue and eigenvector estimates, and quadrature error decays exponentially with node count and minimum eigenvalue distance to the contour.

  • Takeaways & Limitations

    The integral method computes eigenvalues inside a contour without initial eigenvalue or eigenvector approximations and extends to k > m and rank-defective eigenspaces.

  • Takeaways & Limitations

    The method is limited to moderately large problems with a fast sparse direct solver, and a rank-defective example causes two eigenvalues to be missed.

Abstract

from arXiv · show

We propose a numerical method for computing all eigenvalues (and the corresponding eigenvectors) of a nonlinear holomorphic eigenvalue problem that lie within a given contour in the complex plane. The method uses complex integrals of the resolvent operator, applied to at least $k$ column vectors, where $k$ is the number of eigenvalues inside the contour. The theorem of Keldysh is employed to show that the original nonlinear eigenvalue problem reduces to a linear eigenvalue problem of dimension $k$. No initial approximations of eigenvalues and eigenvectors are needed. The method is particularly suitable for moderately large eigenvalue problems where $k$ is much smaller than the matrix dimension. We also give an extension of the method to the case where $k$ is larger than the matrix dimension. The quadrature errors caused by the trapezoid sum are discussed for the case of analytic closed contours. Using well known techniques it is shown that the error decays exponentially with an exponent given by the product of the number of quadrature points and the minimal distance of the eigenvalues to the contour.

1 Introduction

The paper addresses global computation of nonlinear holomorphic eigenvalues inside a contour using resolvent contour integrals and Keldysh’s theorem. The method reduces the problem to a smaller linear eigenproblem and extends to cases with many eigenvalues or dependent eigenvectors.

  • Motivation: Polynomial truncation methods require suitable convergence behavior and may need power-series computations at many points.Their success depends strongly on the radius of convergence and coefficient-matrix decay.
  • Motivation: Newton-type methods can compute eigenvalue groups and associated subspaces, but remain local methods despite robust initial-value convergence.The paper therefore targets the global problem.
  • Contribution: Contour integrals and Keldysh’s theorem provide the basis for computing all eigenvalues and eigenvectors inside a prescribed closed contour.The theorem supplies the resolvent expansion needed by the integral approach.
  • Contribution: For k eigenvalues inside the contour with k ≤ m, the method reduces the nonlinear problem to a linear eigenproblem of dimension k.The contour integrals are applied to a generally random matrix of k columns.
  • Numerical realization: The method is mainly intended for moderately large problems because quadrature with N nodes requires solving Nk linear systems and therefore benefits from a fast sparse direct solver.This requirement is identified as the main numerical effort and a practical scope boundary.
  • Results and extensions: N ≈25 quadrature nodes usually suffice for good eigenvalue and eigenvector estimates, while quadrature error decays exponentially with node count and eigenvalue distance from the contour.An extension handles k > m and rank defects in eigenspaces, including cases with more eigenvalues than the matrix dimension.

2 Nonlinear eigenvalues and Keldysh’ Theorem

This section defines nonlinear eigenvalues and develops Keldysh’s resolvent representation through generalized eigenvectors and adjoint generalized eigenvectors. The resulting residue formulas connect eigenvalue structure to contour-integral computation.

  • Definitions: A nonlinear eigenvalue is a number λ for which T(λ)v = 0 with nonzero v; its vector v is a right eigenvector.The matrix-valued mapping T is assumed holomorphic on an open domain.
  • Resolvent structure: Every eigenvalue is isolated, and T(z) is meromorphic near it with a finite-order pole in the resolvent expansion.The pole order is uniquely determined and is called the order of the pole.
  • Generalized eigenvectors: Root functions and chains of generalized eigenvectors encode higher-order zeros and multiplicities at an eigenvalue.Canonical systems of generalized eigenvectors exist, with ordered partial multiplicities.
  • Keldysh theorem: Keldysh’s theorem represents the singular part of T(z)^−1 using generalized eigenvectors of T and its adjoint T^H.The theorem also identifies a corresponding canonical system for the adjoint problem.
  • Contour formulas: For compact sets or contours containing finitely many eigenvalues, the resolvent admits a holomorphic remainder plus singular contributions from those eigenvalues.The residue theorem then converts these local representations into contour-integral identities.
  • Simple eigenvalues: When eigenvalues are simple, the contour formulas simplify to terms involving left and right eigenvectors normalized through T′(λn).The normalization is expressed as w_n^H T′(λ_n)v_n = 1.

3 The algorithm for a few eigenvalues

The algorithm uses contour integrals of the resolvent to reduce a nonlinear eigenvalue problem inside Γ to a small linear eigenproblem, with numerical quadrature and rank tests enabling eigenvalue and eigenvector recovery.

  • Problem scope: The method targets problems with k eigenvalues inside Γ and, in high-dimensional settings, expects k ≪ m.The initial setup assumes k ≤ m; the opposite case is handled separately.
  • Integral construction: It computes two contour integrals, A0 and A1, using a chosen matrix ˆV and the resolvent evaluations.The integrals are evaluated by solving linear systems at quadrature points.
  • Reduced problem: An SVD of A0 identifies a rank-k subspace, after which a computable k×k matrix B is formed from A0 and A1.The singular values separate the retained rank from numerical zero singular values.
  • Eigenpair recovery: The reduced matrix B is diagonalizable with exactly the eigenvalues inside Γ, and its eigenvectors yield eigenvectors of the nonlinear problem.For simple eigenvalues, eigenvectors are obtained from the eigenvectors of B through the reduced basis V0.
  • Extensions: The extension handles multiple eigenvalues, rank-deficient eigenspaces, and cases with more eigenvalues than the matrix dimension.The generalized method preserves eigenvalues and partial multiplicities but is computationally more expensive in rank-deficient cases.
  • Quadrature and validation: For analytic contours, equidistant trapezoid nodes approximate the integrals, while residual and interior-contour tests validate computed eigenpairs.The main computational cost is solving Nl linear systems, and eigenvalues near Γ can complicate rank testing.

4 Error analysis and numerical examples

The analysis establishes exponential trapezoid-sum convergence for contour-integral computations and evaluates the integral algorithm across polynomial, boundary-value, and rank-defective examples.

  • Error analysis: Trapezoid-sum errors for analytic periodic integrands decay exponentially with the number of quadrature nodes and the strip’s holomorphy width.The contour parametrization is assumed to extend holomorphically to a periodic strip.
  • Error analysis: For a circular contour, the convergence rate depends on the minimum distance d(T) between eigenvalues and the contour.The estimate also covers the case where the spectrum inside the relevant domain is empty.
  • Error analysis: The integral algorithm detects the correct rank for sufficiently large N, enabling perturbation-based eigenvalue error estimates.These conclusions follow from singular-value decomposition and simple-eigenvalue perturbation theory.
  • Example 4.9: Eight eigenvalues inside the test circle are detected and well approximated, with characteristic errors decaying exponentially as N increases.The two tracked eigenvalues exhibit approximately the same decay rate.
  • Example 4.9: Singular-value separation occurs near N ≈25, while the adaptive algorithm reduces retained singular values from 23 at N = 20 to 8 for N ≥95.The fixed-column experiment uses l = 11 columns, whereas the adaptive computation starts with l = 23.
  • Examples 4.10–4.11: For random complex quadratic problems, six interior eigenvalues are well approximated, but nearby exterior eigenvalues weaken singular-value separation without necessarily creating spurious eigenvalues.Very small singular values, approximately 10^-10, can produce spurious eigenvalues and are filtered by the residual test.
  • Example 4.11: In a finite-element nonlinear boundary problem, five eigenvalues in [2, 298] are computed and eigenvector residuals decay exponentially with N.The residual decay is less smooth than in the earlier examples, and the contour avoids a singularity at z = 1.
  • Example 4.11: A rank-defective quadratic problem causes the integral method to miss two of three interior eigenvalues, while only one singular value remains of order one as N increases.The missed eigenvalues share the eigenvector e1, violating the method’s rank condition.

5 The algorithm for many eigenvalues

The extended integral method handles contours containing more eigenvalues than the matrix dimension and rank-deficient eigenspaces by reducing the problem to a k-dimensional linear eigenproblem under rank conditions. Numerical examples show correct eigenvalue recovery, while singular-value and residual tests guide reliability.

  • 5 The algorithm for many eigenvalues: The method is extended to problems with m < k and to rank-deficient eigenspaces.These cases are not covered directly by the basic construction.
  • 5.1 Construction of algorithm: The generalized construction uses additional contour-integral moments, an SVD rank test, and a reduced matrix D whose eigenproblem yields the computed eigenvalues and eigenvectors.The algorithm increases l or K when the numerical rank is insufficient, then accepts eigenpairs using an interior-contour and residual test.
  • 5.1 Construction of algorithm: Under the stated rank conditions, D has the same eigenvalues and partial multiplicities as the nonlinear problem inside Γ.The theorem applies when Γ contains no eigenvalues on its contour and the relevant rank conditions hold.
  • 5.1 Construction of algorithm: The generalization is needed when too many eigenvalues lie inside the contour or when rank defects occur, unlike ordinary linear eigenvalue problems.The paper compares this situation with polynomial linearization, while restricting the extra construction to these nonlinear cases.
  • 5.2 Numerical Examples: In Example 5.4, the algorithm correctly reproduces eigenvalues a = −0.2 and b = 1, with three surviving singular values for K = 2 and l = 3.The example uses a quadratic matrix polynomial with rank defect and compares the results with polyeig.
  • 5.2 Numerical Examples: In Example 5.5, the method finds five eigenvalues inside the delay-system contour, matching previously computed values; smaller N than 150 provides sufficient accuracy.The experiment uses l = 2, K = 3, and N = 150 for the reported computation.
Loading 1003.1580v1…