Measuring the error

All iterative algorithms need some form of termination criteria. In NEP-PACK, all NEP-solvers provide the possibility to specify the desired tolerance, as well as how the error is measured or estimated. The tolerance is specified in the kwarg tol (which is a real number) and the way to measure the error is given in errmeasure.

Standard usage

NEP-PACK comes with several ways to measure errors for many NEP-types.

errmeasure=ResidualErrmeasure(nep): The error is estimated by the use of the residual norm:

\[\mathrm{err}=\frac{\|M(λ)v\|}{\|v\|}.\]

errmeasure=StandardSPMFErrmeasure(nep): The error is estimated by using backward error theory. This error measure will not work for all NEPs. An implementation is provided for any AbstractSPMF. If your NEP is an AbstractSPMF with terms:
\[M(λ)=A_1f_1(λ)+\cdots+A_mf_m(λ)\]
the error will be estimated by
\[\mathrm{err}=\frac{\|M(λ)v\|}{\|v\|}\frac{1}{\|A_1\|_F|f_1(λ)|+\cdots+\|A_m\|_F|f_m(λ)|}.\]
In other words, the StandardSPMFErrmeasure is a weighting of the ResidualErrmeasure.
errmeasure=DefaultErrmeasure(nep): When this errmeasure is specified, NEP-PACK tries to determine a error measure for you. In general, StandardSPMFErrmeasure will be preferred if possible. This behavior may change in future versions of NEP-PACK.
errmeasure=EigvalReferenceErrmeasure(nep,λref): This errmeasure is used when an exact (or very accurate) eigenvalue is already known. Typically, if you wish to visualize the eigenvalue error of a specific method, you run the method twice and use the result of the first run as to instantiate this error measure and get real eigenvalue errors as output.
errmeasure=(λ,v)-> compute_error(λ,v): A user defined error measure can be specified using a function. The function should be take an eigenpair as input, and return a real value. See ErrmeasureType for an example.

Example: Most NEP-solvers take the errmeasure as an kwarg.

julia> nep=nep_gallery("qdep0");
julia> # Solve the problem to residual norm 1e-8
julia> (λ,v)=mslp(nep,errmeasure=ResidualErrmeasure(nep),tol=1e-8);
julia> norm(compute_Mlincomb(nep,λ,v))/norm(v) # It's smaller than tol?
3.503700738108789e-9
julia> nep isa AbstractSPMF # Is it an AbstractSPMF so we can use StandardSPMFErrmeasure?
true
julia> (λ,v)=mslp(nep,errmeasure=StandardSPMFErrmeasure(nep),tol=1e-10);
julia> fv = get_fv(nep); Av = get_Av(nep);
julia> factor=abs(fv[1](λ))*norm(Av[1])+
     abs(fv[2](λ))*norm(Av[2])+abs(fv[3](λ))*norm(Av[3]);
julia> norm(compute_Mlincomb(nep,λ,v))/(norm(v)*factor)
1.6591695084414612e-11

User defined error measure

A common situation is that you want to report the error (as a function of iteration) with a reference solution. We take this situation as an example and show how a relative error of the eigenvalue estimate, compared to a reference solution, can be computed. Compare with type EigvalReferenceErrmeasure which measures an absolute error relative a reference solution.

There are two ways that a user can specify how to measure the error.

User defined error 1: Function handle

The user can provide a function handle which is called to obtain the error. The errmeasure can be a function, which takes two parameters as input (λ,v) and returns the error (or estimate thereof).

If we want to get a very accurate approximation of the true error, we can run the algorithm twice, and the second time we run the algorithm we use the result of the first run as a reference.

julia> nep=nep_gallery("qdep0");
julia> v0=ones(size(nep,1));
julia> (λref,_)=resinv(nep,v=v0,λ=-0.1,logger=0);
julia> myerrmeasure = (λ,v) -> abs(λ-λref)/abs(λ);
julia> (λ,v)=resinv(nep,v=v0,λ=-0.1,logger=1,tol=1e-10,errmeasure=myerrmeasure);
Precomputing linsolver
iter 1 err:0.02854168838549373 λ=-0.1 + 0.0im
iter 2 err:0.8397508140476416 λ=-0.6418389474323298 + 0.0im
iter 3 err:0.17336372619725743 λ=-0.08765753239354723 + 0.0im
iter 4 err:0.0005771170619943501 λ=-0.1029135620110966 + 0.0im
iter 5 err:4.762006833879597e-7 λ=-0.10285411985934721 + 0.0im
iter 6 err:4.074039107701665e-7 λ=-0.10285421074175707 + 0.0im
iter 7 err:2.6448037288912206e-8 λ=-0.10285417155884034 + 0.0im
iter 8 err:1.3926542408883378e-9 λ=-0.10285416898178967 + 0.0im
iter 9 err:6.324560618281378e-11 λ=-0.10285416884505445 + 0.0im

User defined error 2: A user defined type

Due to the multiple dispatch and handling of types in Julia, code may run faster if one uses types instead of function handles. It is possible to do the same simulation as above with a user defined type.

You first need to define a new type

julia> struct RefErrmeasure <: Errmeasure; end

The error measure should then provided in the function estimate_error which we now define as the relative eigenvalue error:

julia> function NonlinearEigenproblems.estimate_error(e::RefErrmeasure,λ,v)
         return abs(λ-λref)/abs(λ);
       end
julia> (λ,v)=resinv(nep,v=v0,λ=0.1,logger=1,tol=1e-10,errmeasure=RefErrmeasure());
Precomputing linsolver
iter 1 err:0.02854168838549373 λ=-0.1 + 0.0im
iter 2 err:0.8397508140476416 λ=-0.6418389474323298 + 0.0im
iter 3 err:0.17336372619725743 λ=-0.08765753239354723 + 0.0im
iter 4 err:0.0005771170619943501 λ=-0.1029135620110966 + 0.0im
iter 5 err:4.762006833879597e-7 λ=-0.10285411985934721 + 0.0im
iter 6 err:4.074039107701665e-7 λ=-0.10285421074175707 + 0.0im
iter 7 err:2.6448037288912206e-8 λ=-0.10285417155884034 + 0.0im
iter 8 err:1.3926542408883378e-9 λ=-0.10285416898178967 + 0.0im
iter 9 err:6.324560618281378e-11 λ=-0.10285416884505445 + 0.0im

As a NEP-solver developer

NEP-solvers should use the Errmeasure as follows. The NEP-solver should take as input an object of the type Errmeasure or function. The fact that it can be different types, is transparent and a NEP-solver developer does not have to do anything to take care of that if the following procedure is followed.

Suppose your solver is defined in a function with this signature:

function mysolver(nep::NEP; errmeasure::ErrmeasureType=DefaultErrmeasure(nep), tol::Real=eps()*100)

In the main for loop you want to call the estimate_error function:

for k=1:maxit
    err=estimate_error(errmeasure,λ,v)
    if (err < tol)
       return (λ,v)
    end
    ....

end

Methods and types

NonlinearEigenproblems.NEPTypes.Errmeasure — Type.

abstract type Errmeasure; end

Concrete subtypes of Errmeasure represent specific ways of measuring the error of an eigenpair. NEP-solvers take such an object as input. As a NEP-solver user, you use the type as follows

julia> quasinewton(nep,errmeasure=ResidualErrmeasure(nep))

User-specified ways of measuring error can be given by creating a new subtype of Errmeasure and using it as a errmeasure keyword. You need to specify the way to measure the error in the method estimate_error.

Note that in practice a Function can essentially be used instead of a Errmeasure-object, which is a simple way to have user-defined error measures. See ErrmeasureType.

NonlinearEigenproblems.NEPTypes.DefaultErrmeasure — Type.

struct DefaultErrmeasure <: Errmeasure
function DefaultErrmeasure(nep::NEP)

When you specify this Errmeasure, NEP-PACK tries to determine a suitable Errmeasure based on the type of the NEP. Note that this behavior may change in future versions.

See also: Errmeasure

NonlinearEigenproblems.NEPTypes.StandardSPMFErrmeasure — Type.

struct StandardSPMFErrmeasure <: Errmeasure
function StandardSPMFErrmeasure(nep::AbstractSPMF)

This Errmeasure provides a way to compute the backward error. The backward error estimate are only given for NEPs which are subtypes of AbstractSPMF. We use the Frobenius norm as the matrix norm, since it is much cheaper to compute than the spectral norm.

Example

julia> nep=nep_gallery("qdep0");
julia> (λ,v)=quasinewton(nep,λ=-1,v=ones(size(nep,1)),errmeasure=StandardSPMFErrmeasure(nep),tol=1e-10,logger=1);
Precomputing linsolver
iter 1 err:0.022010375110869937 λ=-1.0 + 0.0im
iter 2 err:0.002515422247048546 λ=-0.7063330111559607 + 0.0im
iter 3 err:0.000892354247568813 λ=-0.8919579082730457 + 0.0im
iter 4 err:5.445678793151584e-5 λ=-1.0097584042560848 + 0.0im
iter 5 err:6.649967517409105e-7 λ=-1.0023823873044 + 0.0im
iter 6 err:1.0557281809769784e-8 λ=-1.0024660870524031 + 0.0im
iter 7 err:6.420125566431444e-9 λ=-1.0024677891861997 + 0.0im
iter 8 err:3.181093707909799e-10 λ=-1.0024669496893164 + 0.0im
iter 9 err:2.6368050026394416e-11 λ=-1.0024669918249076 + 0.0im

See also: Errmeasure

NonlinearEigenproblems.NEPTypes.ResidualErrmeasure — Type.

struct ResidualErrmeasure <: Errmeasure
function ResidualErrmeasure(nep::NEP)

This Errmeasure species that the residual norm should be used to measure the error.

See also: Errmeasure

NonlinearEigenproblems.NEPTypes.EigvalReferenceErrmeasure — Type.

struct EigvalReferenceErrmeasure{X<:Number} <: Errmeasure
function EigvalReferenceErrmeasure(nep,λref)

Use the difference between a precomputed λ-value (reference solution) and the eigenvalue estimate as the error measure.

Example

julia> using LinearAlgebra
julia> nep=nep_gallery("qdep0");
julia> (λref,vref)=quasinewton(nep,λ=-1,v=ones(size(nep,1)));
julia> (λ,v)=quasinewton(nep,errmeasure=EigvalReferenceErrmeasure(nep,λref),λ=-1.0 ,logger=1,tol=5e-13,v=ones(size(nep,1)));
Precomputing linsolver
iter 1 err:0.002466988585763996 λ=-1.0 + 0.0im
iter 2 err:0.2961339774298033 λ=-0.7063330111559607 + 0.0im
iter 3 err:0.11050908031271833 λ=-0.8919579082730457 + 0.0im
iter 4 err:0.007291415670320767 λ=-1.0097584042560848 + 0.0im
iter 5 err:8.460128136400513e-5 λ=-1.0023823873044 + 0.0im
iter 6 err:9.015333608530796e-7 λ=-1.0024660870524031 + 0.0im
iter 7 err:8.006004357241636e-7 λ=-1.0024677891861997 + 0.0im
iter 8 err:3.889644761834177e-8 λ=-1.0024669496893164 + 0.0im
iter 9 err:3.2391436199930013e-9 λ=-1.0024669918249076 + 0.0im
iter 10 err:2.418474309706653e-10 λ=-1.0024669883439166 + 0.0im
iter 11 err:2.0230705999324528e-11 λ=-1.0024669886059947 + 0.0im
iter 12 err:0.0 λ=-1.002466988585764 + 0.0im

See also: Errmeasure

NonlinearEigenproblems.NEPTypes.estimate_error — Function.

function estimate_error(E::ErrmeasureType,λ,v)

Returns the error estimate for the eigenpair (λ,v). The way to measure the error is specified in E, which can be an Errmeasure or a Function.

NonlinearEigenproblems.NEPTypes.ErrmeasureType — Constant.

ErrmeasureType = Union{Type{<:Errmeasure}, Function}

The ErrmeasureType represents (essentially) what you can insert in the errmeasure keyword argument for most NEP-solvers. It can be a function or an Errmeasure object. If it is a Function this function will be used to obtain error estimate.

Example

This shows how to compute a reference solution and then use this as a reference solution. The error in the second run will be effectively the eigenvector error (appropriately normalized).

julia> using LinearAlgebra
julia> nep=nep_gallery("qdep0");
julia> (λref,vref)=quasinewton(nep,λ=-1,v=ones(size(nep,1)));
julia> myerrmeasure=(λ,v) -> norm(vref/vref[1]-v/v[1]);
julia> (λ,v)=quasinewton(nep,errmeasure=myerrmeasure,λ=-1.0 ,logger=1,tol=5e-13,v=ones(size(nep,1)));
Precomputing linsolver
iter 1 err:46.40296482739195 λ=-1.0 + 0.0im
iter 2 err:2.1592671533657883 λ=-0.7063330111559607 + 0.0im
iter 3 err:0.17079231439255405 λ=-0.8919579082730457 + 0.0im
iter 4 err:0.1633846991066227 λ=-1.0097584042560848 + 0.0im
iter 5 err:0.003434042059262583 λ=-1.0023823873044 + 0.0im
iter 6 err:0.0003182517281689052 λ=-1.0024660870524031 + 0.0im
iter 7 err:2.0105257231740345e-5 λ=-1.0024677891861997 + 0.0im
iter 8 err:1.618661190265619e-6 λ=-1.0024669496893164 + 0.0im
iter 9 err:1.233489068442819e-7 λ=-1.0024669918249076 + 0.0im
iter 10 err:9.44707811957546e-9 λ=-1.0024669883439166 + 0.0im
iter 11 err:7.867601351698812e-10 λ=-1.0024669886059947 + 0.0im
iter 12 err:0.0 λ=-1.002466988585764 + 0.0im

The eigenvalue error can be measured with the EigvalReferenceErrmeasure.

See also: Errmeasure