Projection

struct NewtonInnerSolver
function NewtonInnerSolver(;tol=1e-13,maxit=80,starting_vector=:Vk,
                           newton_function=augnewton)

Uses a Newton-like method to solve the inner problem, with tolerance, and maximum number of iterations given by tol and maxit. The starting vector can be :ones, :randn, or :Vk. The value :Vk specifies the use of the outer NEP-solver keyword argument (Vk). This is typically the previous iterate in the outer method.

The kwarg newton_function, specifies a Function which is called. The type supports augnewton, newton, resinv quasinewton, newtonqr. In principle it can be any method which takes the same keyword arguments as these methods.

See also: InnerSolver, inner_solve

struct IARInnerSolver
function IARInnerSolver(;tol=1e-13,maxit=80,
           starting_vector=:ones,normalize_DEPs=:auto,
           iar_function=iar)

Uses iar, tiar or iar_chebyshev to solve the inner problem, with tolerance, and maximum number of iterations given by tol and maxit. The starting vector can be :ones or :randn. The iar_function can be iar, tiar or iar_chebyshev (or any function taking the same parameters as input). normalize_DEPs determines if the we should carry out precomputation and make sure the projection of a DEP is again a DEP (which can speed up performance). It can take the value true, false or :auto. :auto sets it to true if we use the iar_chebyshev solver.

The kwarg iar_function, specifies a Function which is called. Examples of functions are iar and iar_chebyshev. It can be any NEP-solver which takes the same keyword arguments as these methods.

See also: InnerSolver, inner_solve

function IARChebInnerSolver(;tol=1e-13,maxit=80,starting_vector=:ones,
                            normalize_DEPs=true)

Uses iar_chebyshev to solve the inner problem. See IARInnerSolver for keyword argument documentation.

See also: IARInnerSolver, InnerSolver, inner_solve

struct ContourBeynInnerSolver <: InnerSolver
function ContourBeynInnerSolver(;tol=sqrt(eps(real(Float64))),
                                radius=:auto,N=1000)

Uses contour_beyn to solve the inner problem, with radius and number of quadrature nodes, given by radius and n. If the variable radius is set to :auto, the integration radius will be automatically selected by using the eigenvalue approximations specified by the outer solver.

See also: InnerSolver, inner_solve

struct PolyeigInnerSolver <: InnerSolver
function PolyeigInnerSolver()

Specifies the use of polyeig to solve the inner problem. This is intended for NEPs of the type PEP.

See also: InnerSolver, inner_solve

struct SGIterInnerSolver <: InnerSolver

Uses sgiter to solve the inner problem.

See also: InnerSolver, inner_solve

struct NleigsInnerSolver <: InnerSolver
function NleigsInnerSolver(;Σ= :auto,nodes =:auto, tol=1e-6 )

Uses nleigs to solve the inner problem, in the region Σ with shifts nodes and with tolerance tol. If the variable Σ is set to :auto, the region Σ will be set by using the eigenvalues approximations. See nleigs for description of parameters.

See also: InnerSolver, inner_solve

struct DefaultInnerSolver <: InnerSolver

Dispatches a version of inner_solve based on the type of the NEP provided. This function tries to automatically detect which solver is recommended.

See also: InnerSolver, inner_solve

abstract type InnerSolver

Structs inheriting from this type are used to solve inner problems in an inner-outer iteration.

The InnerSolver objects are passed to the NEP-algorithms, which uses it to dispatch the correct version of the function inner_solve. Utilizes existing implementations of NEP-solvers and inner_solve acts as a wrapper to these.

Example

There is a DefaultInnerSolver that dispatches an inner solver based on the provided NEP. However, this example shows how you can force nlar to use the IARInnerSolver for a PEP.

julia> nep=nep_gallery("pep0", 100);
julia> λ,v = nlar(nep, inner_solver_method=NEPSolver.IARInnerSolver(), neigs=1, num_restart_ritz_vecs=1, maxit=70, tol=1e-8);
julia> norm(compute_Mlincomb(nep,λ[1],vec(v)))
8.68118417430353e-9

See also: inner_solve, DefaultInnerSolver, NewtonInnerSolver, PolyeigInnerSolver, IARInnerSolver, IARChebInnerSolver, SGIterInnerSolver, ContourBeynInnerSolver

inner_solve(is::InnerSolver,T_arit,nep;kwargs...)

Solves the projected linear problem with solver specied with is. This is to be used as an inner solver in an inner-outer iteration. T specifies which method to use. The most common choice is DefaultInnerSolver. The function returns (λv,V) where λv is an array of eigenvalues and V a matrix with corresponding vectors. The struct T_arit defines the arithmetics used in the outer iteration and should prefereably also be used in the inner iteration.

Different inner_solve methods take different kwargs. These are standardized kwargs:

• neigs: Number of wanted eigenvalues (but less or more may be returned)
• σ: target specifying where eigenvalues
• λv, V: Vector/matrix of guesses to be used as starting values
• j: the jth eigenvalue in a min-max characterization
• tol: Termination tolarance for inner solver
• inner_logger: Determines how the inner solves are logged. See Logger for further references

pnep=create_proj_NEP(orgnep::ProjectableNEP[,maxsize [,T]])

Create a NEP representing a projected problem $N(λ)=W^HM(λ)V$, where the base NEP is represented by orgnep. The optional parameter maxsize::Int determines how large the projected problem can be and T is the Number type used for the projection matrices (default ComplexF64). These are needed for memory preallocation reasons. Use set_projectmatrices! and expand_projectmatrices! to specify projection matrices $V$ and $W$.

Example:

The following example illustrates that a projection of a NEP is also a NEP and we can for instance call compute_Mder on it:

julia> nep=nep_gallery("pep0")
julia> V=Matrix(1.0*I,size(nep,1),2);
julia> W=Matrix(1.0*I,size(nep,1),2);
julia> pnep=create_proj_NEP(nep);
julia> set_projectmatrices!(pnep,W,V);
julia> compute_Mder(pnep,3.0)
2×2 Array{Complex{Float64},2}:
  6.08082+0.0im  -5.47481+0.0im
 0.986559+0.0im  -6.98165+0.0im
julia> W'*compute_Mder(nep,3.0)*V  # Gives the same result
2×2 Array{Float64,2}:
 6.08082   -5.47481
 0.986559  -6.98165

If you know that you will only use real projection matrices, you can specify this in at the creation:

julia> pnep=create_proj_NEP(nep,2,Float64);
julia> set_projectmatrices!(pnep,W,V);
julia> compute_Mder(pnep,3.0)
2×2 Array{Float64,2}:
 6.08082   -5.47481
 0.986559  -6.98165

set_projectmatrices!(pnep::Proj_NEP,W,V)

Set the projection matrices for the NEP to W and V, i.e., corresponding the NEP: $N(λ)=W^HM(λ)V$. See also create_proj_NEP.

Example

This illustrates if W and V are vectors of ones, the projected problem becomes the sum of the rows and columns of the original NEP.

julia> nep=nep_gallery("pep0")
julia> pnep=create_proj_NEP(nep);
julia> V=ones(200,1);  W=ones(200,1);
julia> set_projectmatrices!(pnep,W,V);
julia> compute_Mder(pnep,0)
1×1 Array{Complex{Float64},2}:
 48.948104019482756 + 0.0im
julia> sum(compute_Mder(nep,0),dims=[1,2])
1×1 Array{Float64,2}:
 48.948104019482955

expand_projectmatrices!(nep::Proj_SPMF_NEP, Wnew, Vnew)

The projected NEP is updated by adding the last column of Wnew and Vnew to the basis. Note that Wnew and Vnew contain also the "old" basis vectors. See also create_proj_NEP

Example:

In the following example you see that the expanded projected problem has one row and column more, and the leading subblock is the same as the smaller projected NEP.

julia> nep=nep_gallery("pep0"); n=size(nep,1);
julia> V=Matrix(1.0*I,n,2); W=Matrix(1.0*I,n,2);
julia> pnep=create_proj_NEP(nep);
julia> set_projectmatrices!(pnep,W,V);
julia> compute_Mder(pnep,0)
2×2 Array{Complex{Float64},2}:
 0.679107+0.0im   -0.50376+0.0im
 0.828413+0.0im  0.0646768+0.0im
julia> Vnew=[V ones(n)]
julia> Wnew=[W ones(n)]
julia> expand_projectmatrices!(pnep,Wnew,Vnew);
julia> compute_Mder(pnep,0)
3×3 Array{Complex{Float64},2}:
 0.679107+0.0im   -0.50376+0.0im  -12.1418+0.0im
 0.828413+0.0im  0.0646768+0.0im   16.3126+0.0im
 -17.1619+0.0im   -10.1628+0.0im   48.9481+0.0im

abstract ProjectableNEP <: NEP

A ProjectableNEP is a NEP which can be projected, i.e., one can construct the problem $W'*M(λ)Vw=0$ with the Proj_NEP. A NEP which is of this must have the function create_proj_NEP(orgnep::ProjectableNEP) implemented. This function must return a Proj_NEP.

See also: set_projectmatrices!.

Example:

julia> nep=nep_gallery("dep0");
julia> typeof(nep)
DEP{Float64,Array{Float64,2}}
julia> isa(nep,ProjectableNEP)
true
julia> projnep=create_proj_NEP(nep);
julia> e1 = Matrix(1.0*I,size(nep,1),1);
julia> set_projectmatrices!(projnep,e1,e1);
julia> compute_Mder(nep,3.0)[1,1]
-2.942777908030041
julia> compute_Mder(projnep,3.0)
1×1 Array{Complex{Float64},2}:
 -2.942777908030041 + 0.0im

abstract type Proj_NEP <: NEP

Proj_NEP represents a projected NEP. The projection is defined as the NEP

where $M(λ)$ is a base NEP and W and V rectangular matrices representing a basis of the projection spaces. Instances are created with create_proj_NEP. See create_proj_NEP for examples.

Any Proj_NEP needs to implement two functions to manipulate the projection:

• set_projectmatrices!: Set matrices W and V
• expand_projectmatrices!: Effectively expand the matrices W and V with one column.

struct Proj_SPMF_NEP <: Proj_NEP

This type represents the (generic) way to project NEPs which are AbstractSPMF. See examples in create_proj_NEP.

compute_rf(eltype::Type,nep::NEP,x,inner_solver::InnerSolver;
           y=x, target=zero(T), λ0=target,TOL=eps(real(T))*1e3,max_iter=10)

Computes the Rayleigh functional of the nep, i.e., computes a vector $Λ$ of values $λ$ such that $y^TM(λ)x=0$, using the procedure specified in inner_solver. The default behaviour consists of a scalar valued Newton-iteration, and the returned vector has only one element.

The given eltype<:Number is the type of the returned vector.

Example

This uses iar to solve the (scalar) nonlinear problem.

julia> nep=nep_gallery("dep0");
julia> x=ones(size(nep,1));
julia> s=compute_rf(ComplexF64,nep,x,IARInnerSolver())[1] # Take just first element
-1.186623627962043 - 1.5085094961223182im
julia> x'*compute_Mlincomb(nep,s,x)
-8.881784197001252e-16 + 1.0880185641326534e-14im