The NLEVP module

Module containing routines to solve and perturb nonlinear eigenvalue problems.

LinearOperatorFamily

Type for a linear operator family. The type is mutable.

Fields

• terms: Array of terms forming the operator family
• params: Dictionary mapping parameter symbols to their values
• eigval: the symbol of the parameter that is the eigenvalue
• auxval: the symbol of the parameter that is the auxiliary eigenvalue
• active: list of symbols that can be actively change when using an instance of

the family as a function.

• mode: mode that is controlling which terms are evalauted when calling an

instance of the family. This field is not to be modified by the user.

See also: Solution, Term

L=LinearOperatorFamily(fname::String)

Load and construct LinearOperatorFamily from file fname.

See also: save

Solution

Type for the solution of an iterative eigensolver.

Fields

• params: Dictionary mapping parameter symbols to their values
• v: right (direct) eigenvector
• v_adj: left (adjoint) eigenvector
• eigval: the symbol of the parameter that is the eigenvalue
• eigval_pert: extra data for asymptotic series expansion of the eigenvalue
• v_pert: extra data for asymptotic series expansion of the eigenvector
• auxval: the symbol of the parameter that is the auxiliary eigenvalue

See also: LinearOperatorFamily

sol=Solution(params,v,v_adj,eigval,auxval=Symbol())

Standard constructor for a solution object. See Solution for details.

Term{T}

Single term of a linear operator family.

Fields

• coeff::T: matrix coefficient of the term
• func::Tuple: tuple of scalar-valued functions multiplying the matrix coefficient
• symbol::String: character string for displaying the function(s)
• params::Tuple: tuple of tuples containing function symbols for each function
• operator:String: String for displaying the matrix coefficient
• varlist::Array{Symbol,1}: Array containing all symbols of fucntion arguments

that are used at least once

term=Term(coeff,func::Tuple,params::Tuple,operator::String)

Standard constructor for type Term.

Arguments

• coeff: matrix coefficient of the term
• func::Tuple: tuple of scalar-valued functions multiplying the matrix coefficient
• params::Tuple: tuple of tuples containing function symbols for each function
• operator:String: String for displaying the matrix coefficient

Notes

The rendering of the functions of term ist automotized. If the passed functions implement a method f(z::Symbol) then this method is used, otherwise the functions will be simply displayed with the string f. You can overwrite this behavior by explicitly passing a string for the function display using the syntax

term =Term(coeff,func::Tuple,params::Tuple,symbol::String,operator::String)

where the extra argument symbolis the string used to display the function.

See also: LinearOperatorFamily

Ω, P = beyn(L::LinearOperatorFamily, Γ; <keyword arguments>)

Compute all eigenvalues of L inside the contour Γ together with the associated eigenvectors. The contour is given as a list of complex numbers which are interpreted as polygon vertices in the complex plane. The eigenvalues are stored in the list Ω and the eigenvectors are stored in the columns of P. The eigenvector P[:,i] corresponds to the eigenvalue Ω[i], i.e., they satisfy L(Ω[i])P[:,i]=0.

Arguments

• L::LinearOperatorFamily: Definition of the non-linear eigenvalue problem
• Γ::Array: List of complex points defining the contour.
• l::Integer = 5: estimate of the number of eigenvalues inside of Γ.
• K::Integer = 1: Augmention dimension if Γ is assumed to contain more than size(L)[1] eigenvalues.
• N::Integer = 16: Number of evaluation points used to perform the Gauss-Legendre integration along each edge of Γ.
• tol::Float=0.0: Threshold value to discard spurious singular values. If set to 0 (the default) no singular values are discarded.
• pos_test::bool=true: If set to true perform positions test on the computed eigenvalues, i.e., check whether the eigenvalues are enclosed by Γ and disregard all eigenvalues which fail the test.
• output::bool=true: Show progressbar if true.
• random::bool=false: Randomly initialize the V matrix if true.

Returns

• Ω::Array: List of computed eigenvalues
• P::Matrix: Matrix of eigenvectors.

Notes

The original algorithm was first presented by Beyn in [1]. The implementation closely follows the pseudocode from Buschmann et al. in [2].

References

[1] W.-J. Beyn, An integral method for solving nonlinear eigenvalue problems, Linear Algebra and its Applications, 2012, 436(10), p.3839-3863, https://doi.org/10.1016/j.laa.2011.03.030

[2] P.E. Buschmann, G.A. Mensah, J.P. Moeck, Solution of Thermoacoustic Eigenvalue Problems with a Non-Iterative Method, J. Eng. Gas Turbines Power, Mar 2020, 142(3): 031022 (11 pages) https://doi.org/10.1115/1.4045076

See also: inveriter, lancaster, mslp, rf2s, traceiter

A=compute_moment_matrices(L,Γ, <kwargs>)

Compute moment matrices of L integrating along Γ.

Arguments

• L::LinearOperatorFamily
• Γ: List of complex points defining the contour.
• l::Int=5: (optional) estimate of the number of eigenvalues inside of Γ.
• K::Int=1: (optional) Augmention dimension if Γ is assumed to contain more than size(L)[1] eigenvalues.
• N::Int=16: (optional) Number of evaluation points used to perform the Gauss-Legendre integration along each edge of Γ.
• output::Bool=false (optional) toggle waitbar
• random::Bool=false (optional) toggle random initialization of V-matrix

Returns

• A::Array: Moment matrices.

Notes

A is a 3 dimensional array. The last index loops over the various Moments, i.e., A[:,:,i]=∫_Γ z^i inv(L(z)) V dz. If random=false V is initialized with ones on its main diagonal. Otherwise it is random.

n=count_poles_and_zeros(L,Γ;N=16,output=false)

Count number "n" of poles and zeros of the determinant of L inside the contour Γ. The optional parameters N and output. respectively determine the order of the Gauss-Legendre integration and toggle whether output is desplayed.

Notes

Pole orders count negative while zero order count positive. The evaluation is based on applying the residue theorem to the determinant. It utilizes Jacobi's formula to compute the necessary derivatives from a trace operation. and is therefore only suitable for small problems.

msg=decode_error_flag(flag::Int)

Decode error flag flag from an iterative solver into a string msg.

See also: inveriter, lancaster, mslp, rf2s, traceiter

sol::Solution, n, flag = householder(L::LinearOperatorFamily, z; <keyword arguments>)

Use a Householder method to iteratively find an eigenpar of L, starting the the iteration from z.

Arguments

• L::LinearOperatorFamily: Definition of the nonlinear eigenvalue problem.
• z: Initial guess for the eigenvalue.
• maxiter::Integer=10: Maximum number of iterations.
• tol=0: Absolute tolerance to trigger the stopping of the iteration. If the difference of two consecutive iterates is abs(z0-z1)<tol the iteration is aborted.
• relax=1: relaxation parameter
• lam_tol=Inf: tolerance for the auxiliary eigenvalue to test convergence. The default is infinity, so there is effectively no test.
• order::Integer=1: Order of the Householder method, Maximum is 5
• nev::Integer=1: Number of Eigenvalues to be searched for in intermediate ARPACK calls.
• v0::Vector: Initial vector for Krylov subspace generation in ARPACK calls. If not provided the vector is initialized with ones.
• v0_adj::Vector: Initial vector for Krylov subspace generation in ARPACK calls. If not provided the vector is initialized with v0.
• output::Bool: Toggle printing online information.

Returns

• sol::Solution
• n::Integer: Number of perforemed iterations
• flag::Integer: flag reporting the success of the method. Most important values are 1: method converged, 0: convergence might be slow, -1:maximum number of iteration has been reached. For other error codes see the source code.

Notes

Housholder methods are a generalization of Newton's method. If order=1 the Housholder method is identical to Newton's method. The solver then reduces to the "generalized Rayleigh Quotient iteration" presented in [1]. If no relaxation is used (relax == 1), the convergence rate is of order+1. With relaxation (relax != 1) the convergence rate is 1. Thus, with a higher order less iterations will be necessary. However, the computational time must not necessarily improve nor do the convergence properties. Anyway, if the method converges, the error in the eigenvalue is bounded above by tol. For more details on the solver, see the thesis [2].

References

[1] P. Lancaster, A Generalised Rayleigh Quotient Iteration for Lambda-Matrices,Arch. Rational Mech Anal., 1961, 8, p. 309-322, https://doi.org/10.1007/BF00277446

[2] G.A. Mensah, Efficient Computation of Thermoacoustic Modes, Ph.D. Thesis, TU Berlin, 2019

See also: beyn

sol,n,flag = inveriter(L,z;maxiter=10,tol=0., relax=1., x0=[],v=[], output=true)

Deploy inverse iteration for finding an eigenpair of L.

Arguments

• L::LinearOperatorFamily: Definition of the nonlinear eigenvalue problem.
• z: Initial guess for the eigenvalue.
• maxiter::Integer=10: Maximum number of iterations.
• tol=0: Absolute tolerance to trigger the stopping of the iteration. If the difference of two consecutive iterates is abs(z0-z1)<tol the iteration is aborted.
• relax=1: relaxation parameter
• x0::Vector: initial guess for eigenvector. If not provided it is all ones.
• v::Vector: normalization vector. If not provided it is all ones.
• output::Bool: Toggle printing online information.

Returns

• sol::Solution
• n::Integer: Number of perforemed iterations
• flag::Integer: flag reporting the success of the method. Most important values are 1: method converged, 0: convergence might be slow, -1:maximum number of iteration has been reached. For other error codes see the source code.

Notes

The algorithm uses Newton-Raphson itreations on the vector level for iteratively solving the problem L(λ)x == 0 and v'x == 0. The implementation is based on Algorithm 1 in [1]

References

[1] V. Mehrmann and H. Voss (2004), Nonlinear eigenvalue problems: A challenge for modern eigenvalue methods, GAMM-Mitt. 27, 121–152.

See also: beyn, lancaster, mslp, rf2s, traceiter, decode_error_flag

sol,n,flag = lancaster(L,z;maxiter=10,tol=0., relax=1., x0=[],y0=[], output=true)

Deploy Lancaster's Rayleigh-quotient iteration for finding an eigenvalue of L.

Arguments

• L::LinearOperatorFamily: Definition of the nonlinear eigenvalue problem.
• z: Initial guess for the eigenvalue.
• maxiter::Integer=10: Maximum number of iterations.
• tol=0: Absolute tolerance to trigger the stopping of the iteration. If the difference of two consecutive iterates is abs(z0-z1)<tol the iteration is aborted.
• relax=1: relaxation parameter
• x0::Vector: initial guess for right iteration vector. If not provided it is all ones.
• y0::Vector: initial guess for left iteration vector. If not provided it is all ones.
• output::Bool: Toggle printing online information.

Returns

• sol::Solution
• n::Integer: Number of perforemed iterations
• flag::Integer: flag reporting the success of the method. Most important values are 1: method converged, 0: convergence might be slow, -1:maximum number of iteration has been reached. For other error codes see the source code.

Notes

The algorithm is a generalization of Rayleigh-quotient-iteration by Lancaster [1].

References

[1] P. Lancaster, A Generalised Rayleigh Quotient Iteration for Lambda-Matrices,Arch. Rational Mech Anal., 1961, 8, p. 309-322, https://doi.org/10.1007/BF00277446

See also: beyn, inveriter, mslp, rf2s, traceiter, decode_error_flag

Ω,P = moments2eigs(A; tol_σ=0.0, return_σ=false)

Compute eigenvalues Ω and eigenvectors P from moment matricesA`.

Arguments

• A::Array: moment matrices
• tol::Float=0.0: tolerance for truncating singular values.
• return_σ::Bool=false: if true also return the singular values as third return value.

Returns

• Ω::Array: List of computed eigenvalues
• P::Matrix: Matrix of eigenvectors.
• Σ::Array: List of singular values (only if return_σ==true)

Notes

The the i-th coloumn in P corresponds to the i-th eigenvalue in Ω, i.e., Ω[i] and P[:,i] are an eigenpair.

See also: compute_moment_matrices

sol,n,flag = mslp(L,z;maxiter=10,tol=0.,relax=1.,lam_tol=Inf,order=1,nev=1,v0=[],v0_adj=[],num_order=1,scale=1.0+0.0im,output=false)

Deploy the method of successive linear problems for finding an eigentriple of L.

Arguments

• L::LinearOperatorFamily: Definition of the nonlinear eigenvalue problem.
• z: Initial guess for the eigenvalue.
• maxiter::Integer=10: Maximum number of iterations.
• tol=0: Absolute tolerance to trigger the stopping of the iteration. If the difference of two consecutive iterates is abs(z0-z1)<tol the iteration is aborted.
• relax=1: relaxation parameter
• lam_tol=Inf: tolerance for the auxiliary eigenvalue to test convergence. The default is infinity, so there is effectively no test.
• order::Integer=1: Order of the Householder method, Maximum is 5
• nev::Integer=1: Number of Eigenvalues to be searched for in intermediate ARPACK calls.
• v0::Vector: Initial vector for Krylov subspace generation in ARPACK calls. If not provided the vector is initialized with ones.
• v0_adj::Vector: Initial vector for Krylov subspace generation in ARPACK calls. If not provided the vector is initialized with v0.
• num_order::Int=1: order of the numerator polynomial when forming the Padé approximant.
• scale::Number=1: scaling factor applied to the tolerance (and the eigenvalue output when output==true).
• output::Bool=false: Toggle printing online information.

Returns

• sol::Solution
• n::Integer: Number of perforemed iterations
• flag::Integer: flag reporting the success of the method. Most important values are 1: method converged, 0: convergence might be slow, -1:maximum number of iteration has been reached. For other error codes see the source code.

Notes

This is a variation of the method of succesive linear Problems [1] as published in Algorithm 3 in [2]. It treats the eigenvalue ω as a parameter in a linear eigenvalue problem L(ω)p=λYp and computes the root of the implicit relation λ(ω)==0.

The default values "order==1" and "num_order==1" will result in a Newton-iteration for finding the roots. Choosing a higher order will result in Housholder methods as a generalization of Newton's method. If no relaxation is used (relax == 1), the convergence rate is of order+1. With relaxation (relax != 1) the convergence rate is 1. Thus, with a higher order less iterations will be necessary. However, the computational time must not necessarily improve nor do the convergence properties. Anyway, if the method converges, the error in the eigenvalue is bounded above by tol. For more details on the solver, see the thesis [2].

Householder methods are based on expanding the relation λ=λ(ω) into a [1/order-1]-Padé approximant. The numerator order may be increased using the keyword "num_order". This feature is experimental.

References

[1] A.Ruhe, Algorithms for the non-linear eigenvalue problem, SIAM J. Numer. Anal. ,10:674–689, 1973.

[2] V. Mehrmann and H. Voss (2004), Nonlinear eigenvalue problems: A challenge for modern eigenvalue methods, GAMM-Mitt. 27, 121–152.

[3] G.A. Mensah, Efficient Computation of Thermoacoustic Modes, Ph.D. Thesis, TU Berlin, 2019 doi:10.14279/depositonce-8952

See also: beyn, inveriter, lancaster, rf2s, traceiter, decode_error_flag

perturb!(sol::Solution,L::LinearOperatorFamily,param::Symbol,N::Int; <keyword arguments>)

Compute the Nth order power series perturbation coefficients for the solution sol of the nonlineaer eigenvalue problem given by the operator family L with respect to the parameter param. The coefficients will be stored in the field sol.eigval_pert and sol.v_pert for the eigenvalue and the eignevector, respectively.

Keyword Arguments

• mode = :compact: parameter controlling internal programm flow. Use the default, unless you know what you are doing.

Notes

For large perturbation orders N the method might be slow.

See also: perturb_fast!

perturb_fast!(sol::Solution,L::LinearOperatorFamily,param::Symbol,N::Int; <keyword arguments>)

Compute the Nth order power series perturbation coefficients for the solution sol of the nonlineaer eigenvalue problem given by the operator family L with respect to the parameter param. The coefficients will be stored in the field sol.eigval_pert and sol.v_pert for the eigenvalue and the eignevector, respectively.

Keyword Arguments

• mode = :compact: parameter controlling internal programm flow. Use the default, unless you know what you are doing.

Notes

This method reads multi-indeces for the computation of the power series coefficients from disk. Make sure that JulHoltz is properly installed to use this fast method.

See also: perturb!

perturb_norm!(sol::Solution,L::LinearOperatorFamily,param::Symbol,N::Int; <keyword arguments>)

Compute the Nth order power series perturbation coefficients for the solution sol of the nonlineaer eigenvalue problem given by the operator family L with respect to the parameter param. The coefficients will be stored in the field sol.eigval_pert and sol.v_pert for the eigenvalue and the eignevector, respectively.

Keyword Arguments

• mode = :compact: parameter controlling internal programm flow. Use the default, unless you know what you are doing.

Notes

This method reads multi-indeces for the computation of the power series coefficients from disk. Make sure that JulHoltz is properly installed to use this fast method.

See also: perturb!

Ω,P=pos_test(Ω,P,Γ)

Filter eigenpairs Ω and Pfor those whose eigenvalues lie inside Γ.

See also: moments2eigs

sol,n,flag = rf2s(L,z;maxiter=10,tol=0., relax=1., x0=[],y0=[], output=true)

Deploy two sided Rayleigh-functional iteration for finding an eigentriple of L.

Arguments

• L::LinearOperatorFamily: Definition of the nonlinear eigenvalue problem.
• z: Initial guess for the eigenvalue.
• maxiter::Integer=10: Maximum number of iterations.
• tol=0: Absolute tolerance to trigger the stopping of the iteration. If the difference of two consecutive iterates is abs(z0-z1)<tol the iteration is aborted.
• relax=1: relaxation parameter
• x0::Vector: initial guess for right eigenvector. If not provided it is the first basis vector [1 0 0 0 ...].
• y0::Vector: initial guess for left eigenvector. If not provided it is the first basis vector `[1 0 0 0 ...].
• output::Bool: Toggle printing online information.

Returns

• sol::Solution
• n::Integer: Number of perforemed iterations
• flag::Integer: flag reporting the success of the method. Most important values are 1: method converged, 0: convergence might be slow, -1:maximum number of iteration has been reached. For other error codes see the source code.

Notes

The algorithm is known to have a cubic convergence rate. Its implementation follows Algorithm 4.9 in [1].

References

[1] S. Güttel and F. Tisseur, The Nonlinear Eigenvalue Problem, 2017, http://eprints.ma.man.ac.uk/2531/.

See also: beyn, inveriter, lancaster, mslp, traceiter, decode_error_flag

save(fname::String,L::LinearOperatorFamily)

Save L to file fname. The file is utf8 encoded and adheres to a Julia-enriched TOML standard.

See also: LinearOperatorFamily

solve(L,Γ;<keywordarguments>)

Find eigenvalues of linear operator family L inside contour Γ.

Arguments

• L::LinearOperatorFamily: Linear operator family for which eigenvalues are to be found
• Γ::List: List of complex points defining the contour.

Keywordarguments (optional)

• Δl::Int=1: Increment to the search space dimension in Beyn's algorithm
• N::Int=16: Number of integration points per edge of Γ.
• tol::Float=1E-8: Tolerance for local eigenvalue solver.
• eigvals::Dict=Dict(): Dictionary of known eigenvalues and corresponding eigenvectors.
• maxcycles::Int=1: Number of local correction cycles to global eigenvalue estimate.
• nev::Int=1: Number of tested eigenvalues for finding the best local match to the global estimates.
• max_outer_cycles::Int=1: Number of outer cycles, i.e., incremental increases to the search space dimension by Δl.
• atol_σ::Float=1E-12: absolute tolerance for terminating the inner cycle if the maximum singular value is σ<atol_σ`
• rtol_σ::Float=1E-8: relative tolerance for terminating the inner cycle if the the ratio of the initial and the current maximum singular value is σ/σ0<rtol_σ
• loglevel::Int=0: loglevel value between 0 (no output) and 2 (max output) to control the detail of printed output.

Returns

• eigvals::Dict:List of computed eigenvalues and corresponding eigenvectors.

Notes

The algorithm is experimental. It startes with a low-dimensional Beyn integral. From which estimates to the eigenvalues are computed. These estimates are used to analytically correct Beyn's integral. This step requires no additional numerical integration and is therefore relatively quick. New local estimates are computed and if the local solver then found new eigenvalues these are used to further correct Beyn's integral. In an optional outer loop the entire procedure can be repeated, after increasing the search space dimension.

See also: householder, beyn

sol,n,flag = traceiter(L,z;maxiter=10,tol=0.,relax=1.,output=true)

Deploy trace iteration for finding an eigevalue of L.

Arguments

• L::LinearOperatorFamily: Definition of the nonlinear eigenvalue problem.
• z: Initial guess for the eigenvalue.
• maxiter::Integer=10: Maximum number of iterations.
• tol=0: Absolute tolerance to trigger the stopping of the iteration. If the difference of two consecutive iterates is abs(z0-z1)<tol the iteration is aborted.
• relax=1: relaxation parameter
• output::Bool: Toggle printing online information.

Returns

• sol::Solution
• n::Integer: Number of perforemed iterations
• flag::Integer: flag reporting the success of the method. Most important values are 1: method converged, 0: convergence might be slow, -1:maximum number of iteration has been reached. For other error codes see the source code.

Notes

The algorithm applies Newton-Raphson iteration for finding the root of the determinant. The updates for the iterates are computed by using trace operations and Jacobi's formula [1]. The trace operation renders the method slow and inaccurate for medium and large problems.

References

[1] S. Güttel and F. Tisseur, The Nonlinear Eigenvalue Problem, 2017, http://eprints.ma.man.ac.uk/2531/.

See also: beyn, inveriter, lancaster, mslp, rf2s, decode_error_flag

w=wn(z,Γ)

Compute winding number, i.e., how often Γ is wrapped around z. The sign of the winding number indocates the wrapping direction.

d=size(L::LinearOperatorFamily)

Get dimension d of d×d linear operator family L.

V=initialize_V(d::Int,l::Int,random::Bool=false)

initialioze matrix with drows and l colums either as diagonal matrix featuring ones on the main diagonal and anywhere else or with random entries (random=true).

See also: generate_subspace

f=polyval(p,z)

Evaluate polynomial f(z)=∑_i p[i]z^1 at z, using Horner's scheme.

P::LinearOperatorFamily=project(L::LinearOperatorFamily,Q)

Project L on the subspace spanned by the unitary matrix Q, i.e. P(z)=Q'*L(z)*Q.

See also: generate_subspace

This is a minimal TOML parser. See https://en.wikipedia.org/wiki/TOML