NEP-Solvers

The NEP solver methods implemented in NEP-PACK, are accessed by the functions below. The functions all return $λ,v$ where $λ$ is either a number (eigenvalue) a vector of eigenvalues $v$ is either a vector containing an eigenvector or a matrix whose columns corresponding to the eigenvectors. Two-sided methods may return $λ,v,w$ where $w$ are the left eigenvectors.

The first optional parameter in all NEP solver methods is a type. This type specifies which arithmetic should be used for the algorithm.

Example:

julia> nep=nep_gallery("dep0");
julia> λ,v=augnewton(ComplexF64,nep,v=ones(5))
(-0.15955391823299256 + 0.0im, Complex{Float64}[0.12505315954062152 + 0.0im, 0.8475907515488971 + 0.0im, -0.10910413290558324 + 0.0im, 0.027714719799174125 + 0.0im, 0.10874550201689052 + 0.0im])
julia> typeof(λ)
Complex{Float64}
julia> λ,v=augnewton(Float16,nep,v=ones(5))
(Float16(-0.718), Float16[0.435, 0.6606, -0.205, -0.1445, 0.254])
julia> typeof(λ)
Float16

The NEP-solvers can be separated into the following types (with some overlap):

Newton type methods
Projection methods
Contour integral methods
Arnoldi and Krylov based methods
Class specific methods

Newton type methods

NonlinearEigenproblems.NEPSolver.newton — Function.

λ,v = newton([eltype],nep::NEP;[errmeasure,][tol,][maxit,][λ,][v,][c,][logger,][armijo_factor=1,][armijo_max])

Applies Newton-Raphsons method on the system of nonlinear equations with n+1 unknowns:

\[M(λ)v=0\]

\[c^Hv-1=0\]

The vector c is the orthogonalization vector. If c=0 the current approximation will be used for the orthogonalization. See augnewton for other parameters.

Example

julia> using LinearAlgebra
julia> nep=nep_gallery("dep0");
julia> λ,v=newton(nep);
julia> minimum(svdvals(compute_Mder(nep,λ)))
1.9997125567227177e-16

References

Nichtlineare Behandlung von Eigenwertaufgaben, Z. Angew. Math. Mech. 30 (1950) 281-282.
A. Ruhe, Algorithms for the nonlinear eigenvalue problem, SIAM J. Numer. Anal. 10 (1973) 674-689

NonlinearEigenproblems.NEPSolver.augnewton — Function.

augnewton([eltype], nep::NEP; [errmeasure,][tol,][maxit,][λ,][v,][c,][logger,][linsolvercreator,][armijo_factor,][armijo_max])

Run the augmented Newton method. The method is equivalent to newton() in exact arithmetic, but works only with operations on vectors of length n.

The following keyword arguments are in common for many NEP-solvers:

logger is either a Logger object or an Int. If it is an Int, a PrintLogger(logger) will be instantiated. logger=0 prints nothing, logger=1 prints more, etc.
errmeasure determines how error is measured. It is either a function handle or an object of the type Errmeasure. If it is a function handle, it should take (λ,v) as input and return a real scalar (the error). See Errmeasure and ErrmeasureType for further description.
tol is a scalar which determines termination. If errmeasure is less than tol the eigenpair is marked as converged.
The scalar λ and the vector v are starting approximations.
maxit determines the maximum number of iterations. The error NoConvergenceException is thrown if this is exceeded.
The linsolvecreator specifies how the linear system should be solved. See LinSolver for further information.
armijo_factor specifies if an Armijo rule should be applied, and its value specifies the scaling factor of the step length (per reduction step). The variable armijo_max specifies the maximum number of step length reductions.

Example

This illustrates the equivalence between newton and augnewton.

julia> nep=nep_gallery("dep1")
julia> λ1,v1=newton(nep,maxit=20,v=ones(size(nep,1)),λ=0)
julia> λ2,v2=augnewton(nep,maxit=20,v=ones(size(nep,1)),λ=0)
julia> λ1-λ2
0.0 + 0.0im

References

Nichtlineare Behandlung von Eigenwertaufgaben, Z. Angew. Math. Mech. 30 (1950) 281-282.
A. Ruhe, Algorithms for the nonlinear eigenvalue problem, SIAM J. Numer. Anal. 10 (1973) 674-689

NonlinearEigenproblems.NEPSolver.resinv — Function.

λ,v = resinv([eltype],nep::NEP;[errmeasure,][tol,][maxit,][λ,][v,][c,][logger,][armijo_factor=1,][armijo_max,][linsolvecreator])

Applies residual inverse iteration method for nonlinear eigenvalue problems. The kwarg linsolvecreator is a function which specifies how the linear system is created. The function calls compute_rf for the computation of the Rayleigh functional. See augnewton for other parameters.

Example

The example shows how to specify if the method should run in real or complex mode (or any other Number type).

julia> nep=nep_gallery("qdep0");
julia> λ,v=resinv(nep,λ=-2,v=ones(size(nep,1)))
julia> typeof(λ)
Complex{Float64}
julia> norm(compute_Mlincomb(nep,λ,v))
6.688224435370382e-12
julia> λ,v=resinv(Float64,nep,λ=-2,v=ones(size(nep,1)))
julia> typeof(λ)
Float64
julia> norm(compute_Mlincomb(nep,λ,v))
5.939894690000396e-12

References

A. Neumaier, Residual inverse iteration for the nonlinear eigenvalue problem, SIAM J. Numer. Anal. 22 (1985) 914-923

NonlinearEigenproblems.NEPSolver.quasinewton — Function.

quasinewton([T=ComplexF64],nep,[errmeasure,][tol,][maxit,][λ,][v][ws][logger][linsolvercreator,][armijo_factor,][armijo_max])

An implementation of the quasi-Newton approach referred to as quasi-Newton 2 in the reference. The method involves one linear system solve per iteration corresponding with the matrix $M(λ)$, where $λ$ is constant. The vector ws is a representation of the normalization, in the sense that $c^T=w_s^TM(λ)$, where all iterates satisfy $c^Tx_i=1$. See augnewton for other parameters.

Example

julia> nep=nep_gallery("pep0")
julia> λ,v=quasinewton(nep,λ=1.0,v=ones(size(nep,1)));
julia> norm(compute_Mlincomb(nep,λ,v))/norm(v)
5.448264607410413e-12

References

Jarlebring, Koskela, Mele, Disguised and new Quasi-Newton methods for nonlinear eigenvalue problems, Numer. Algorithms, 79:311-335, 2018. preprint

NonlinearEigenproblems.NEPSolver.mslp — Function.

 mslp([eltype],nep::NEP;[errmeasure,][tol,][maxit,][λ,][logger,][eigsolvertype::Type])

Runs the method of successive linear problems. The method requires the solution of a generalized eigenvalue problem in every iteration. The method used for the eigenvalue computation is specified in eigsolvertype. See augnewton for other parameters.

Example

Create a rational NEP using a SPMF_NEP.

julia> eye=Matrix{Float64}(I,3,3);
julia> Av=[ones(3,3),eye,triu(ones(3,3))];
julia> fv=[S-> S, S -> S^2, S->inv(S-one(S)*10)];
julia> nep=SPMF_NEP(Av,fv);
julia> (λ,v)=mslp(nep);
julia> compute_Mlincomb(nep,λ,v)
3-element Array{Complex{Float64},1}:
 -1.38778e-17+1.65715e-18im
 -5.55112e-17+1.30633e-17im
 -4.16334e-17-1.54436e-17im

References

A. Ruhe, Algorithms for the nonlinear eigenvalue problem, SIAM J. Numer. Anal. 10 (1973) 674-689

NonlinearEigenproblems.NEPSolver.sgiter — Function.

λ,v = sgiter([eltype],nep::NEP,j::Integer;[λ_min,][λ_max,][λ,][errmeasure,][tol,][maxit,][logger,][eigsolvertype::Type,])

Finds the j-th eigenvalue of the NEP using safeguarded iteration, with eigenvalue numbering according to min-max theory. The method only works for Hermitian problems, and the eigenvalues are assumed to be real. If an interval [λ_min,λ_max] is given, then the Rayleigh functional is assumed to be unique on the interval. If no interval is given, then the minimum solution is always taken. The method requires the computation of (all) eigenvalues of a matrix. The eigsolvertype is a Type that specifies which eigevalue solver is used inside the algorithm.

See augnewton for other parameters.

Example

julia> nep = nep_gallery("real_quadratic");
julia> λ,v = sgiter(nep, 1, λ_min = -10, λ_max = 0,  λ = -10, maxit = 100);
julia> minimum(svdvals(compute_Mder(nep,λ)))
0.0
julia> norm(v)
1.0

References

V. Mehrmann and H. Voss, Nonlinear eigenvalue problems: a challenge for modern eigenvalue methods, GAMM‐Mitteilungen 27.2 (2004): 121-152.
H. Voss and B. Werner, Solving sparse nonlinear eigenvalue problems. Technical Report 82/4, Inst. f. Angew. Mathematik, Universität Hamburg, 1982.
B. Werner. Das Spektrum von Operatorenscharen mit verallgemeinerten Rayleighquotienten. PhD thesis, Fachbereich Mathematik, Universität Hamburg, 1970

NonlinearEigenproblems.NEPSolver.rfi — Function.

rfi(nep,nept,[λ=0,][errmeasure,][tol=eps()*100,][maxit=100,][v=randn,][u=randn,][logger=0,][linsolvecreator=DefaultLinSolverCreator(),])

This is an implementation of the two-sided Rayleigh functional Iteration (RFI) to compute an eigentriplet of the problem specified by nep. This method requires the transpose of the NEP, specified in nept. λ, u and v are initial guesses for the eigenvalue, the right eigenvector and the left eigenvector respectively. See augnewton for other parameters.

Example

julia> nep=nep_gallery("dep0");
julia> nept=DEP([nep.A[1]',nep.A[2]'])
julia> λ,v,u=rfi(nep,nept)
julia> compute_resnorm(nep,λ,v) # v is a right eigenvector
3.0171586304599647e-16
julia> compute_resnorm(nept,λ,u) # u is a right eigenvector
7.145081514857076e-17

Reference

Algorithm 4 in Schreiber, Nonlinear Eigenvalue Problems: Newton-type Methods and Nonlinear Rayleigh Functionals, PhD thesis, TU Berlin, 2008.

NonlinearEigenproblems.NEPSolver.rfi_b — Function.

rfi_b(nep,nept,[λ=0,][errmeasure,][tol=eps()*100,][maxit=100,][v=randn,][u=randn,][logger=0,][linsolvecreator=DefaultLinSolverCreator(),])

This is an implementation of the two-sided Rayleigh functional Iteration(RFI)-Bordered version to compute an eigentriplet of the problem specified by nep. This method requires the transpose of the NEP, specified in nept. λ, u and v are initial guesses for the eigenvalue, the right eigenvector and the left eigenvector respectively. See augnewton for other parameters.

Example

julia> nep=nep_gallery("dep0");
julia> nept=DEP([nep.A[1]',nep.A[2]'])
julia> λ,v,u=rfi_b(nep,nept,v=ones(ComplexF64,size(nep,1)))
julia> compute_resnorm(nep,λ,v) # v is a right eigenvector
5.343670589284583e-15

Reference

Algorithm 5 in Schreiber, Nonlinear Eigenvalue Problems: Newton-type Methods and Nonlinear Rayleigh Functionals, PhD thesis, TU Berlin, 2008.

NonlinearEigenproblems.NEPSolver.blocknewton — Function.

(S,X)=blocknewton(nep [S,] [X,] [errmeasure,] [tol,] [maxit,] [armijo_factor,] [armijo_max,] [logger])

Applies the block Newton method to nep::AbstractSPMF. The method computes an invariant pair (S,X) using the block Newton approach of Kressner. The variables S,X correspond to starting approximations. The function errmeasure shoule be defined for errmeasure(S,X) and meausures the error in the pair (S,X). See augnewton for other parameters.

Example

The example shows that compute_MM() becomes zero when a solution has been computed.

julia> nep=nep_gallery("dep0",3);
julia> (S,X)= blocknewton(nep)
julia> compute_MM(nep,S,X)
3×2 Array{Complex{Float64},2}:
 -1.11022e-16+0.0im  1.11022e-16+0.0im
          0.0+0.0im          0.0+0.0im
  1.38778e-17+0.0im  2.77556e-17+0.0im

This example solves the gun problem from the Berlin-Manchester collection

julia> using NonlinearEigenproblems.Gallery
julia> nep=nep_gallery("nlevp_native_gun");
julia> II=[1.0 0; 0 1]; S=150^2*II; V=[II;zeros(size(nep,1)-2,2)];
julia> (Z,X)=blocknewton(nep,S=S,X=V,logger=1,armijo_factor=0.5,maxit=20)
Iteration 1: Error: 6.081316e+03
Iteration 2: Error: 1.701970e-02 Armijo scaling=0.031250
Iteration 3: Error: 1.814887e-02 Armijo scaling=0.250000
...
Iteration 13: Error: 6.257442e-09
Iteration 14: Error: 2.525942e-15

References

D. Kressner A block Newton method for nonlinear eigenvalue problems, Numer. Math., 114 (2) (2009), pp. 355-372

NonlinearEigenproblems.NEPSolver.newtonqr — Function.

λ,v = newtonqr([eltype],nep::NEP;[errmeasure,][tol,][maxit,][λ,][v,][c,][logger])

This function implements the Newton-QR method as formulated in the reference. The method ivolves the computation of a rank-revealing QR factorization of $M(λ)$, with the idea that on convergence the the last diagonal element $R[n,n]$ of the upper-triangular matrix $R$ becomes zero as a result of $M(λ)$ becoming singular. Since the computation of a QR factorization is expensive, it is advisable to use this method for problems of small size or problems with a certain structure that makes the QR computation less expensive. See augnewton for other parameters.

Example

julia> nep=nep_gallery("pep0")
julia> λ,v=newtonqr(nep,v=ones(size(nep,1)));
julia> norm(compute_Mlincomb(nep,λ,v))/norm(v)
8.440206093655014e-15

References

Kublanovskaya, V. N., (1970). On an approach to the solution of the generalized latent value problem for λ-matrices, SIAM J. Numer. Anal. 7, 532–537
Güttel, S., & Tisseur, F. (2017). The nonlinear eigenvalue problem. Acta Numerica, 26, 1-94. doi:10.1017/S0962492917000034

NonlinearEigenproblems.NEPSolver.implicitdet — Function.

λ,v = implicitdet([eltype],nep::NEP;[errmeasure,][tol,][maxit,][λ,][v,][c,][logger])

This function implements the Implicit determinant method as formulated Algorithm 4.3 in the reference. The method applies Newton-Raphson to the equation $det(M(λ))/det(G(λ)) = 0$, where $G(λ)$ is a saddle point matrix with $M(λ)$ in the (1,1) block. The (2,1) and (1,2) blocks of $G(λ)$ are set to $c^H$ and $c$ respectively. Note that $G(λ)$ can be non-singular even when $M(λ)$ is singular. See reference for more information. See augnewton for other parameters.

Example

julia> nep=nep_gallery("pep0")
julia> λ,v=implicitdet(nep,v=ones(size(nep,1)));
julia> norm(compute_Mlincomb(nep,λ,v))/norm(v)
2.566371972986362e-14

References

Spence, A., & Poulton, C. (2005). Photonic band structure calculations using nonlinear eigenvalue techniques, J. Comput. Phys., 204 (2005), pp. 65–8
Güttel, S., & Tisseur, F. (2017). The nonlinear eigenvalue problem. Acta Numerica, 26, 1-94. doi:10.1017/S0962492917000034

NonlinearEigenproblems.NEPSolver.broyden — Function.

S,V = broyden([eltype,]nep::NEP[,approxnep];kwargs)

Runs Broydens method (with deflation) for the nonlinear eigenvalue problem defined by nep. An approximate nep can be provided which is used as an initialization of starting matrix/vectors. The optional argument approxnep determines how to initiate the algorithm. It can be an NEP, the symbol :eye corresponding to starting with an identity matrix, and a Matrix (of size $n imes n$). Beside most of the standard kwargs as described in augnewton, it supports pmax which is subspace used in deflation, essentially the number of eigenvalues, add_nans::Bool which determines if NaNs should be added in book keeping. eigmethod which can be :eig, :eigs or :invpow. The :invpow is an implementation of the power method, which is slow but works well e.g. for BigFloat. The kwarg recompute_U determines if the U-matrix should be recomputed in every deflation (which can be more robust). The implementation has two loggers logger and inner_logger. The logger corresponds to outer iterations (deflation) and inner_logger is the iterates in Broydens method. The kwarg check_error_every and print_error_every detemine how often errors should be check and how often they should be printed. For real problems with complex conjugate symmetry, you may want to set the kwarg addconj=true in order to reduce computation by automatically adding the complex conjugate vectors.

The method computes an invariant pair and can therefore find several eigenvalues. The retured value is (S,V) is an invariant pair and the eigenvalues are on the diagonal of S. Eigenpairs can be directly extracted with get_deflated_eigpairs.

Example

julia> nep=nep_gallery("dep0");
julia> S,V=broyden(nep);
julia> λ=S[1,1]
-0.15955391823299253 - 3.874865266487398e-19im
julia> minimum(svdvals(compute_Mder(nep,λ)))
1.6293996560844023e-16
julia> λ=S[2,2]
-0.5032087003825461 + 1.1969823800738464im
julia> minimum(svdvals(compute_Mder(nep,λ)))
1.1073470346550144e-15
julia> λ=S[3,3]
1.2699713558173726 + 5.342786996459857e-16im
julia> minimum(svdvals(compute_Mder(nep,λ)))
5.905315846211231e-16
julia> broyden(nep,logger=2,check_error_every=1);  # Prints out a lot more convergence info

In order to extract eigenpairs you can use the following:

julia> (D,X)=get_deflated_eigpairs(S,V,size(nep,1));
julia> for i=1:3; @show norm(compute_Mlincomb(nep,D[i],X[:,i])); end
norm(compute_Mlincomb(nep, D[i], X[:, i])) = 8.459878994614521e-13
norm(compute_Mlincomb(nep, D[i], X[:, i])) = 1.2102336671048442e-13
norm(compute_Mlincomb(nep, D[i], X[:, i])) = 2.1012363973403225e-16

References

Jarlebring, Broyden’s method for nonlinear eigenproblems, SIAM J. Sci. Comput., 41:A989–A1012, 2019, https://arxiv.org/pdf/1802.07322

Projection methods

NonlinearEigenproblems.NEPSolver.nlar — Function.

function nlar([eltype],nep::ProjectableNEP,[orthmethod=ModifiedGramSchmidt],[neigs=10],[errmeasure],[tol=eps(real(T))*100],[maxit=100],[λ0=0],[v0=randn(T,size(nep,1))],[logger=0],[linsolvercreator=DefaultLinSolverCreator()],[R=0.01],[eigval_sorter=residual_eigval_sorter],[qrfact_orth=false],[max_subspace=100],[num_restart_ritz_vecs=8],[inner_solver_method=DefaultInnerSolver(),][inner_logger=0])

The function implements the Nonlinear Arnoldi method, which finds neigs eigenpairs (or throws a NoConvergenceException) by projecting the problem to a subspace that is expanded in the course of the algorithm. The basis is orthogonalized either by using the QR method if qrfact_orth is true or else by an orthogonalization method orthmethod). This entails solving a smaller projected problem using a method specified by inner_solver_method. The logging of the inner solvers are descided by inner_logger, which works in the same way as logger. (λ0,v0) is the initial guess for the eigenpair. linsolvercreator specifies how the linear system is created and solved. R is a parameter used by the function specified by eigval_sorter to reject those ritz values that are within a distance R from any of the converged eigenvalues, so that repeated convergence to the same eigenpair can be avoided. max_subspace is the maximum allowable size of the basis befor the algorithm restarts using a basis made of num_restart_ritz_vecs ritz vectors and the eigenvectors that the algorithm has converged to.

See augnewton for other parameters.

Example

julia> nep=nep_gallery("dep0_tridiag");
julia> λ,v=nlar(nep,tol=1e-7,neigs=1,maxit=100,v=ones(size(nep,1)));
julia> norm(compute_Mlincomb(nep,λ[1],v))
8.00192341259751e-7

References

H. Voss, An Arnoldi method for nonlinear eigenvalue problems. BIT. Numer. Math. 44: 387-401, 2004.

NonlinearEigenproblems.NEPSolver.jd_betcke — Function.

jd_betcke([eltype]], nep::ProjectableNEP; [neigs=1], [tol=eps(real(T))*100], [maxit=100], [λ=zero(T)], [orthmethod=DGKS],  [errmeasure], [linsolvercreator=DefaultLinSolverCreator()], [v = randn(size(nep,1))], [logger=0], [inner_logger=0], [inner_solver_method=DefaultInnerSolver()], [projtype=:PetrovGalerkin], [target=zero(T)])

The function computes eigenvalues using Jacobi-Davidson method, which is a projection method. The projected problems are solved using a solver spcified through the type inner_solver_method. The logging of the inner solvers are descided by inner_logger, which works in the same way as logger. For numerical stability the basis is kept orthogonal, and the method for orthogonalization is specified by orthmethod, see the package IterativeSolvers.jl. The function tries to compute neigs number of eigenvalues, and throws a NoConvergenceException if it cannot. The value λ and the vector v are initial guesses for an eigenpair. linsolvercreator is a function which specifies how the linear system is created and solved. The target is the center around which eiganvlues are computed. By default the method uses a Petrov-Galerkin framework, with a trial (left) and test (right) space, hence $W^H T(λ) V$ is the projection considered. By specifying projtype to be :Galerkin then W=V.

See augnewton for other parameters.

Example

julia> nep=nep_gallery("dep0",50);
julia> λ,v=jd_betcke(nep,tol=1e-8,maxit=20);
julia> norm(compute_Mlincomb(nep,λ[1],v[:,1]))
1.9570222439914163e-10

References

T. Betcke and H. Voss, A Jacobi-Davidson-type projection method for nonlinear eigenvalue problems. Future Gener. Comput. Syst. 20, 3 (2004), pp. 363-372.
H. Voss, A Jacobi–Davidson method for nonlinear eigenproblems. In: International Conference on Computational Science. Springer, Berlin, Heidelberg, 2004. pp. 34-41.

Contour integral methods

NonlinearEigenproblems.NEPSolver.contour_beyn — Function.

λv,V=contour_beyn([eltype,] nep [,mintegrator];[tol,][logger,][σ,][radius,][linsolvercreator,][N,][neigs,][k])

The function computes eigenvalues using Beyn's contour integral approach, using an ellipse centered at σ with radii given in radius, or if only one radius is given, the contour is a circle. The numerical quadrature method is specified in mintegrator, which is a type inheriting from MatrixIntegrator, by default MatrixTrapezoidal. For a parallell implementation of the integrator use MatrixTrapezoidalParallel. The integer k specifies size of the probe subspace. N corresponds to the number of quadrature points. Ellipses are the only supported contours. The linsolvercreator must create a linsolver that can handle (rectangular) matrices as right-hand sides, not only vectors. We integrate in complex arithmetic so eltype must be complex type.

The kwargs neigs specifies the number of wanted eigvals, and k is the number of columns in the matrix to be integrated (default k=neigs+1). If you give the k parameter and set neigs=typemax(Int) all found eigenvalues will be returned. The kwarg sanity_check decides if sorting and checking (and removal) of eigpairs should be done. If disabled, the method returns k (potentially inaccurate) eigpairs. The parameters errmeasure and tol and rank_drop_tol are used for the sanity check, to extract accurate eigenvalues.

Example

julia> using LinearAlgebra
julia> nep=nep_gallery("dep0");
julia> # Look for two eigvals in unit disk
julia> λv,V=contour_beyn(nep,radius=1.1,neigs=3);
julia> norm(compute_Mlincomb(nep,λv[1],V[:,1])) # Eigenpair 1
1.7462847531404259e-15
julia> norm(compute_Mlincomb(nep,λv[2],V[:,2])) # Eigenpair 2
7.69695692032292e-15

References

Wolf-Jürgen Beyn, An integral method for solving nonlinear eigenvalue problems, Linear Algebra and its Applications 436 (2012) 3839–3863

NonlinearEigenproblems.NEPSolver.contour_block_SS — Function.

contour_block_SS([eltype,] nep [,mintegrator];[tol,][logger,][σ,][radius,][linsolvercreator,][N,][neigs,][k,][L])

This is an implementation of the block_SS contour integral method which is based on the computation of higher order moments. The contour is an ellipse centered at σ with radii given in radius, or if only one radius is given, the contour is a circle. The numerical quadrature method is specified in mintegrator, which is a type inheriting from MatrixIntegrator, by default MatrixTrapezoidal. For a parallell implementation of the integrator use MatrixTrapezoidalParallel. The integer k specifies size of the probe subspace. N corresponds to the number of quadrature points. The integer L specifies the number of moments. Ellipses are the only supported contours. The linsolvercreator must create a linsolver that can handle (rectangular) matrices as right-hand sides, not only vectors. We integrate in complex arithmetic so eltype must be complex type.

Example

julia> nep=SPMF_NEP([[0 1 ; 1 1.0], [1 0 ; 0 0]], [s->one(s),s->exp(1im*s^2)]);
julia> λ,V=contour_assu(nep,radius=3,neigs=6)
julia> @show λ
6-element Array{Complex{Float64},1}:
  4.496403249731884e-15 + 2.506628274630998im
        -2.506628274631 - 2.8727020762175925e-15im
  3.219972424519104e-16 - 2.5066282746310034im
     2.5066282746310096 - 1.1438072192922029e-15im
 -2.3814273710772784e-7 - 7.748469160458366e-8im
   2.381427350935646e-7 + 7.748467479992284e-8im

References

Asakura, Sakurai, Tadano, Ikegami, Kimura, A numerical method for nonlinear eigenvalue problems using contour integrals, JSIAM Letters, 2009 Volume 1 Pages 52-55
Van Beeumen, Meerbergen, Michiels. Connections between contour integration and rational Krylov methods for eigenvalue problems, 2016, TW673, https://lirias.kuleuven.be/retrieve/415487/

Arnoldi and Krylov based methods

IAR

The Infinite ARnoldi method.

NonlinearEigenproblems.NEPSolver.iar — Function.

iar(nep,[maxit=30,][σ=0,][γ=1,][linsolvecreator=DefaultLinSolverCreator(),][tol=eps()*10000,][neigs=6,][errmeasure,][v=rand(size(nep,1),1),][logger=0,][check_error_every=1,][orthmethod=DGKS,][proj_solve=false,][inner_solver_method=DefaultInnerSolver(),][inner_logger=0])

Run the infinite Arnoldi method on the nonlinear eigenvalue problem stored in nep.

The target σ is the center around which eiganvalues are computed. The value γ corresponds to scaling and specifying a shift and scaling is effectively the same as the transformation λ=γs+σ where s is now the eigenvalue parameter. If you want eigenvalues in a disk centered, select σ as the center of the disk and γ as the radius. The vector v is the starting vector for constructing the Krylov space. The orthogonalization method, used in contructing the orthogonal basis of the Krylov space, is specified by orthmethod, see the package IterativeSolvers.jl. The iteration is continued until neigs Ritz pairs have converged. This function throws a NoConvergenceException if the wanted eigenpairs are not computed after maxit iterations. However, if neigs is set to Inf the iteration is continued until maxit iterations without an error being thrown. The parameter proj_solve determines if the Ritz paris are extracted using the Hessenberg matrix (false), or as the solution to a projected problem (true). If true, the method is descided by inner_solver_method, and the logging of the inner solvers are descided by inner_logger, which works in the same way as logger.

See augnewton for other parameters.

Example

julia> using NonlinearEigenproblems, LinearAlgebra
julia> nep=nep_gallery("dep0",100);
julia> v0=ones(size(nep,1));
julia> λ,v=iar(nep;v=v0,tol=1e-5,neigs=3);
julia> norm(compute_Mlincomb!(nep,λ[1],v[:,1])) # Is it an eigenvalue?
julia> λ    # print the computed eigenvalues
3-element Array{Complex{Float64},1}:
 -0.15606211475666945 - 0.12273439802763578im
 -0.15606211475666862 + 0.12273439802763489im
  0.23169243065648365 - 9.464790582509696e-17im

References

Algorithm 2 in Jarlebring, Michiels Meerbergen, A linear eigenvalue algorithm for the nonlinear eigenvalue problem, Numer. Math, 2012

IAR Chebyshev

A Chebyshev version of the IAR method.

NonlinearEigenproblems.NEPSolver.iar_chebyshev — Function.

iar_chebyshev(nep,[maxit=30,][σ=0,][γ=1,][linsolvecreator=DefaultLinSolverCreator(),][tolerance=eps()*10000,][neigs=6,][errmeasure,][v=rand(size(nep,1),1),][logger=0,][check_error_every=1,][orthmethod=DGKS][a=-1,][b=1,][compute_y0_method=ComputeY0ChebAuto])

Run the infinite Arnoldi method (Chebyshev version) on the nonlinear eigenvalue problem stored in nep.

The target σ is the center around which eiganvalues are computed. A Ritz pair λ and v is flagged a as converged (to an eigenpair) if errmeasure is less than tol. The vector v is the starting vector for constructing the Krylov space. The orthogonalization method, used in contructing the orthogonal basis of the Krylov space, is specified by orthmethod, see the package IterativeSolvers.jl. The iteration is continued until neigs Ritz pairs converge. This function throws a NoConvergenceException if the wanted eigenpairs are not computed after maxit iterations. However, if neigs is set to Inf the iteration is continued until maxit iterations without an error being thrown. The kwarg compute_y0_method specifying how the next vector of the Krylov space (in Chebyshev format) can be computed. See compute_y0_cheb in the module NEPSolver with the command ?NEPSolver.compute_y0_cheb.

See augnewton for other parameters.

Example

julia> using NonlinearEigenproblems, LinearAlgebra
julia> nep=nep_gallery("dep0",100);
julia> v0=ones(size(nep,1));
julia> λ,v=iar_chebyshev(nep;v=v0,tol=1e-5,neigs=3);
julia> norm(compute_Mlincomb!(nep,λ[1],v[:,1])) # Is it an eigenvalue?
julia> λ    # print the computed eigenvalues
3-element Array{Complex{Float64},1}:
julia> norm(compute_Mlincomb(nep,λ[1],v[:,1]))
 0.050462487848960284 - 1.4289626573515395e-18im
 -0.07708779190301127 + 7.703053374113074e-18im
   0.1503856540695659 - 1.662582577182149e-17im

References

Algorithm 2 in Jarlebring, Michiels Meerbergen, A linear eigenvalue algorithm for the nonlinear eigenvalue problem, Numer. Math, 2012

For the iar_chebyshev the following compute_y0_cheb method is needed, in order to avoid explicit conversions between the Chebyshev basis and the monimial basis.

NonlinearEigenproblems.NEPSolver.compute_y0_cheb — Function.

y0 = compute_y0_cheb([eltype],nep::NEPTypes.DEP,::Type{ComputeY0ChebPEP},X,Y,M0inv,precomp::AbstractPrecomputeData)

Computes the vector y0 used in iar_chebyshev given by

\[ y_0 = \sum_{i=1}^N T_{i-1}(γ) x_i - \sum_{j=1}^m A_j \left( \sum_{i=1}^{N+1} T_{i-1}(-ρ \tau_j+γ) y_i \right )\]

where T(c) is the vector containing $T_i(c)$ as coefficients, where $T_i$ is the i-th Chebyshev polynomial of the first kind.

y0 = compute_y0_cheb([eltype],nep::NEPTypes.PEP,::Type{ComputeY0ChebPEP},X,Y,M0inv,precomp::AbstractPrecomputeData)

Computes the vector y0 used in iar_chebyshev given by

\[ y_0 = \sum_{j=0}^{d-1} A_{j+1} x D^j T(c) - y T(c)\]

where T(c) is the vector containing $T_i(c)$ as coefficients, where $T_i$ is the i-th Chebyshev polynomial of the first kind and $D$ is the derivation matrix in Chebyshev basis.

y0 = compute_y0_cheb([eltype],nep::NEPTypes.SPMF_NEP,::Type{ComputeY0ChebPEP},X,Y,M0inv,precomp::AbstractPrecomputeData)

Computes the vector y0 used in iar_chebyshev given by

\[ y_0= \sum_{j=0}^{m} M^{(j)}(\mu) X b_j \left( D_N \right) T_N(c) - Y T_N(c)\]

where T(c) is the vector containing $T_i(c)$ as coefficients, where $T_i$ is the i-th Chebyshev polynomial of the first kind and $b_j(\lambda)=(f_j(0)-f_j(\lambda))/\lambda=f[\lambda,0]$ are divided differences.

y0 = compute_y0_cheb([eltype],nep::NEPTypes.NEP,::Type{ComputeY0ChebNEP},X,Y,M0inv,precomp::AbstractPrecomputeData)

Computes the vector y0 used in iar_chebyshev defined as

\[ y_0 =\left( \sum_{i=0}^{N-1} B \left( \frac{d}{d \theta} \right) \hat T_i(\theta) x_i \right)(0) - \sum_{i=0}^{N} T_i(c) y_i\]

where $T_i$ is the i-th Chebyshev polynomial of the first kind, $ \ hat T_i$ is the i-th Chebyshev polynomial of the first kind for the interval [a,b]. For a generic nep, this quantity is computed by converting polynomials in monomial basis. This procedure may be numerical unstable if many iterations are required. If for the specific nep a closed formula is available, we suggest to overload this function.

TIAR

The Tensor Infinite ARnoldi method.

NonlinearEigenproblems.NEPSolver.tiar — Function.

tiar(nep,[maxit=30,][σ=0,][γ=1,][linsolvecreator=DefaultLinSolverCreator(),][tolerance=eps()*10000,][neigs=6,][errmeasure,][v=rand(size(nep,1),1),][logger=0,][check_error_every=1,][orthmethod=DGKS,][proj_solve=false,][inner_solver_method=DefaultInnerSolver(),][inner_logger=0])

Run the tensor infinite Arnoldi method on the nonlinear eigenvalue problem stored in nep. This is equivalent to iar, but handles orthogonalization with a tensor representation.

See augnewton for other parameters.

Example

julia> using NonlinearEigenproblems, LinearAlgebra
julia> nep=nep_gallery("dep0",100);
julia> v0=ones(size(nep,1));
julia> λ,v=tiar(nep;v=v0,tol=1e-5,neigs=3);
julia> norm(compute_Mlincomb!(nep,λ[1],v[:,1])) # Is it an eigenvalue?
julia> λ    # print the computed eigenvalues
3-element Array{Complex{Float64},1}:
 0.050462487743188206 - 3.4059600538119376e-18im
 -0.07708769561361105 + 8.611006691570004e-19im
   0.1503916927814904 + 9.388210527944734e-18im

References

Algorithm 2 in Jarlebring, Mele, Runborg, The Waveguide Eigenvalue Problem and the Tensor Infinite Arnoldi Method, SIAM J. Scient. computing, 39 (3), A1062-A1088, 2017

Infinite Lanczos based methods

The Infinite Bi-Lanczos method.

NonlinearEigenproblems.NEPSolver.infbilanczos — Function.

λv,V,U=infbilanczos([eltype],nep, nept,[linsolvecreator,][linsolvertcreator,][v,][u,][σ,][γ,][tol,][neigs,][errmeasure,][logger,][maxit,][check_error_every])

Executes the Infinite Bi-Lanczos method on the problem defined by nep::NEP and nept::NEP. nep:NEP is the original nonlinear eigenvalue problem and nept::NEP is its (hermitian) transpose: $M(λ^*)^H$. v and u are starting vectors, σ is the shift and γ the scaling. The iteration is continued until neigs Ritz pairs have converged. This function throws a NoConvergenceException if the wanted eigenpairs are not computed after maxit iterations. However, if neigs is set to Inf the iteration is continued until maxit iterations without an error being thrown. See augnewton for other parameters.

Example:

julia> nep=nep_gallery("dep0");
julia> A=get_Av(nep); fv=get_fv(nep);
julia> At=[copy(A[1]'),copy(A[2]'),copy(A[3]')]
julia> nept=SPMF_NEP(At,fv); # Create the transposed NEP
julia> λv,V=infbilanczos(nep,nept,neigs=3,v=ones(size(nep,1)))
julia> norm(compute_Mlincomb(nep,λv[1],V[:,1]))
9.907299130783851e-15

References:

The infinite bi-Lanczos method for nonlinear eigenvalue problems, S. W. Gaaf and E. Jarlebring, SIAM J. Sci. Comput. 39:S898-S919, 2017, preprint

The Infinite Lanczos method, for symmetric NEPs

NonlinearEigenproblems.NEPSolver.ilan — Function.

ilan(nep,[maxit=30,][σ=0,][γ=1,][linsolvecreator=DefaultLinSolverCreator(),][tolerance=eps()*10000,][neigs=6,][errmeasure,][v=rand(size(nep,1),1),][logger=0,][check_error_every=30,][orthmethod=DGKS])

Run the infinite Lanczos method on the symmetric nonlinear eigenvalue problem stored in nep. The current implementation supports only neps in SPMF format.

The target σ is the center around which eiganvalues are computed. The kwarg errmeasure is a function handle which can be used to specify how the error is measured to be used in termination (default is absolute residual norm). A Ritz pair λ and v is flagged a as converged (to an eigenpair) if errmeasure is less than tol. The vector v is the starting vector for constructing the Krylov space. The orthogonalization method, used in contructing the orthogonal basis of the Krylov space, is specified by orthmethod, see the package IterativeSolvers.jl. The iteration is continued until neigs Ritz pairs have converged. This function throws a NoConvergenceException if the wanted eigenpairs are not computed after maxit iterations. However, if neigs is set to Inf the iteration is continued until maxit iterations without an error being thrown.

See augnewton for other parameters.

Example

julia> using NonlinearEigenproblems, LinearAlgebra
julia> nep=nep_gallery("dep_symm_double",10);
julia> v0=ones(size(nep,1));
julia> λ,v=ilan(nep;v=v0,tol=1e-5,neigs=3);
julia> norm(compute_Mlincomb!(nep,λ[1],v[:,1])) # Is it an eigenvalue?
julia> λ    # print the computed eigenvalues
3-element Array{Complex{Float64},1}:
  0.03409997385842267 - 1.425205509434608e-19im
 -0.03100798730589012 + 5.66201058098364e-20im
  -0.0367653644764646 - 1.607494907684445e-19im

References

Algorithm 2 in Mele, The infinite Lanczos method for symmetric nonlinear eigenvalue problems, https://arxiv.org/abs/1812.07557, 2018

NLEIGS

NonlinearEigenproblems.NEPSolver.nleigs — Function.

nleigs(nep::NEP, Σ::AbstractVector{Complex{T}})

Find a few eigenvalues and eigenvectors of a nonlinear eigenvalue problem, using the nleigs algorithm.

Arguments

nep: An instance of a nonlinear eigenvalue problem.
Σ: A vector containing the points of a polygonal target set in the complex plane.
Ξ: A vector containing a discretization of the singularity set.
logger: Level of display (0, 1, 2).
maxdgr: Max degree of approximation.
minit: Min number of iterations after linearization is converged.
maxit: Max number of total iterations.
tol: Tolerance for residual.
tollin: Tolerance for convergence of linearization.
v: Starting vector.
errmeasure: Function for error measure (residual norm). Called with arguments (λ,v).
isfunm : Whether to use matrix functions.
static: Whether to use static version of NLEIGS.
leja: Use of Leja-Bagby points (0 = no, 1 = only in expansion phase, 2 = always).
nodes: Prefixed interpolation nodes (only when leja is 0 or 1).
reusefact: Reuse of matrix factorizations (0 = no, 1 = only after converged linearization, 2 = always).
blksize: Block size for pre-allocation.
return_details: Whether to return solution details (see NleigsSolutionDetails).
check_error_every: Check for convergence / termination every this number of iterations.

See augnewton for other parameters.

Return values

λ: Vector of eigenvalues of the nonlinear eigenvalue problem NLEP inside the target set Σ.
X: Corresponding matrix of eigenvectors.
res: Corresponding residuals.
details: Solution details, if requested (see NleigsSolutionDetails).

Example

julia> nep=nep_gallery("dep0");
julia> unit_square = float([1+1im, 1-1im, -1-1im,-1+1im])
julia> (λ,v)=nleigs(nep,unit_square);
julia> norm(compute_Mlincomb(nep,λ[1],v[:,1]))
5.028313023882492e-14
julia> norm(compute_Mlincomb(nep,λ[2],v[:,2]))
1.1937025845487509e-13

References

S. Guettel, R. Van Beeumen, K. Meerbergen, and W. Michiels. NLEIGS: A class of fully rational Krylov methods for nonlinear eigenvalue problems. SIAM J. Sci. Comput., 36(6), A2842-A2864, 2014.
NLEIGS Matlab toolbox

Class specific methods

Companion linearizations

NonlinearEigenproblems.NEPSolver.companion — Function.

E,A = companion(nep::PEP);

Linearizes a polynomial eigenvalue problem (PEP) a to the companion form, as in the paper by Mehrmann and Voss. More precisely, for a k-th degree PEP with n-by-n coefficient matrices, this returns matrices E and A, both kn-by-kn, corresponding to the linearized problem

\[Ax = λEx\]

Example

julia> pep = nep_gallery("pep0");
julia> E,A = companion(pep);
julia> λ, V = eigen(A,E);
julia> minimum(svd(compute_Mder(pep,λ[1])).S)
2.4845217786736996e-12

References

V. Mehrmann and H. Voss, Non-linear eigenvalue problems, a challenge for modern eigenvalue methods, GAMM‐Mitteilungen (2004)

NonlinearEigenproblems.NEPSolver.polyeig — Function.

λ,v = polyeig([eltype],nep::PEP,[eigsolvertype,])

Linearizes a polynomial eigenvalue problem (PEP) a to the companion form and solves the corresponding linear eigenvalue problem; see companion. The eigsolvertype is optinal can be used to specify how the linear problem is solved; see eig_solve, and EigSolver.

Example

julia> pep = nep_gallery("pep0");
julia> λ,V = polyeig(pep);
julia> minimum(svd(compute_Mder(pep,λ[1])).S)
1.2050763381899922e-14
julia> norm(compute_Mlincomb(pep,λ[2],vec(V[:,2])))
1.0245470569036458e-12

λ,v = polyeig([eltype],nep::ChebPEP)

Computes a companion linearization for the NEP represented in a Chebyshev basis, and returns eigenpairs.

Example

julia> using LinearAlgebra
julia> nep=nep_gallery("dep0");
julia> chebpep=ChebPEP(nep,9,-3,1,cosine_formula_cutoff=5);
julia> (λv,V)=polyeig(chebpep);
julia> ii=argmin(abs.(λv));
julia> λ=λv[ii];
julia> v=V[:,ii];
julia> norm(compute_Mlincomb(chebpep,λ,v))
1.3543968603949142e-14
julia> # Actually, it's not a bad approx to the original NEP either
julia> norm(compute_Mlincomb(nep,λ,v))
4.326355966047557e-6