Types & Data structures

type DEP <: AbstractSPMF
function DEP(AA::Vector{AbstractMatrix} [,tauv::Vector=[0,1.0]])

A DEP (Delay Eigenvalue problem) is a problem defined by the sum

where all of the matrices are of size $n×n$. This type of NEP describes the stability of time-delay systems.

The construction takes the system matrices $A_i$, and tauv is a vector of the values $τ_i$.

Example:

julia> A0=randn(3,3); A1=randn(3,3);
julia> tauv=[0,0.2] # Vector with delays
julia> dep=DEP([A0,A1],tauv)
julia> λ=3.0;
julia> M1=compute_Mder(dep,λ)
julia> M2=-λ*I+A0+A1*exp(-tauv[2]*λ)
julia> norm(M1-M2)
0.0

struct PEP <: AbstractSPMF
function PEP(AA::Vector{AbstractMatrix})

The type PEP defines a polynomial eigenvalue problem via its monomial coefficients. A polynomial eigenvalue problem (PEP) is defined by the sum the

where $i = 0,1,2,$, and all of the matrices are of size $n×n$. The vector AA contains $A_1,...$.

Example

julia> A0=[1.0 3; 4 5]; A1=A0.+one(2); A2=ones(2,2);
julia> pep=PEP([A0,A1,A2])
julia> compute_Mder(pep,3)-(A0+A1*3+A2*9)
2×2 Array{Float64,2}:
 0.0  0.0
 0.0  0.0

See also polyeig, companion, ChebPEP, interpolate.

ChebPEP(orgnep::NEP,k,[a=-1,[b=1]] [,cosine_formula_cutoff=5])

The type ChebPEP<:AbstractSPMF represents a polynomial function where the function is stored using a Chebyshev basis scaled to the interval [a,b], i.e.,

where $T_i$ are the scaled and shifted Chebyshev polynomials.

The constructor ChebPEP takes nep::NEP as an input and interpolates this NEP in k Chebyshev nodes, resulting in a polynomial of degree k-1, represented by its coefficients in the Chebyshev basis. Interpolation in Chebyshev nodes and representation with Chebyshev basis, is known to have attractive approximation properties, as well as robustness with respect to round-off errors.

The kwarg cosine_formula_cutoff decides how the Chebyshev polynomials should be computed. For larger degrees, it is better to use the cosine formula, whereas for low degrees the explicit monomial expression is more efficient. The explicit monomial expression will be used for degrees lower than cosine_formula_cutoff.

Example:

julia> nep=nep_gallery("dep0");
julia> chebpep=ChebPEP(nep,9);
julia> using LinearAlgebra;
julia> norm(compute_Mder(nep,0.3)-compute_Mder(chebpep,0.3))
1.2881862971045282e-8
julia> chebpep=ChebPEP(nep,19); # Better interpolation
julia> norm(compute_Mder(nep,0.3)-compute_Mder(chebpep,0.3))
2.0312004517316714e-15

See also: polyeig, PEP

function REP(A,roots,poles)

A REP-call creates a rational eigenvalue problem. The REP is defined by the sum $Σ_i A_i s_i(λ)/q_i(λ)$, where i = 0,1,2,..., all of the matrices are of size $n×n$ and $s_i$ and $q_i$ are polynomials. The constructor takes the roots and poles as input of polynomials with normalized highest coefficient. The NEP is defined as

where p has the roots roots and q has the roots poles.

Example

julia> A0=[1 2; 3 4]; A1=[3 4; 5 6];
julia> nep=REP([A0,A1],[1,3], [4,5,6]);
julia> compute_Mder(nep,3)
2×2 Array{Float64,2}:
 Inf  Inf
 Inf  Inf
julia> (λ,x)=quasinewton(nep,v=[1;0])
(-0.3689603779201249 + 0.0im, Complex{Float64}[-2.51824+0.0im, 1.71283+0.0im])
julia> -λ*x+A0*x+A1*x*(λ-1)*(λ-3)/((λ-4)*(λ-5)*(λ-6))
2-element Array{Complex{Float64},1}:
 -2.5055998942313806e-13 + 0.0im
   1.318944953254686e-13 + 0.0im

abstract NEP

A NEP object represents a nonlinear eigenvalue problem. All NEPs should implement

size(nep::NEP,d)

and at least one of the following

• M = compute_Mder(nep::NEP,λ::Number,i::Integer=0)
• V = compute_Mlincomb(nep::NEP,λ::Number,V::AbstractVecOrMat,a::Vector) (or compute_Mlincomb!)
• MM = compute_MM(nep::NEP,S,V)

struct SPMF_NEP{T<:AbstractMatrix,Ftype}  <: AbstractSPMF{T}
function SPMF_NEP(AA, fii [,check_consistency=true] [,Schur_fact = false]
                  [,align_sparsity_patterns = false] [,Ftype=ComplexF64])

An SPMF_NEP is a NEP defined by a Sum of Products of Matrices and Functions, i.e.,

All of the matrices $A_0,...$ are of size $n×n$ and $f_i$ are a functions. The functions $f_i$ must be defined for matrices in the standard matrix function sense. The constructor creates a SPMF_NEP consisting of matrices AA and functions fii.

Parameters

• AA is a Vector of matrices. The matrices have to be of the same type. If you need a NEP with different types you can use SumNEP to construct a sum of two SPMF_NEP.

• fii is a Vector of functions. Each function takes one parameter S. The functions must be available both as a scalar valid function and a matrix function. If S is a square matrix, fii[k](S) musst also be a square matrix. If S is a scalar fii[k](S) is a scalar.

• check_consistency (default true) determines if we should initiate by running tests to verify that the fii satisfies the conditions that every function is valid both for matrices and scalars. This is done by using @code_typed and the functions need to be type-stable in that sense.

• align_sparsity_patterns (default false) has effect only for sparse matrices (SparseMatrixCSC). If align_sparsity_patterns=true the SparseMatrixCSC matrices will be replaced by equivalent SparseMatrixCSC matrices where the colptr and rowval are identical. This increases the speed of some functions, e.g., compute_Mder. If align_sparsity_patterns=true the matrices in the NEP should be considered read only. If the sparsity patterns are completely or mostly distinct, it may be more efficient to set this flag to false.

• Ftype (default ComplexF64) determines an underlying type of the functions. The output of any function should be "smaller" than the promoted type of the input and Ftype. More precisely, if F=fii[k], then the type logic is as follows eltype(F(λ))=promote_type(eltype(λ),Ftype).

• Schur_fact (default false) determines if the compute_MM function should triangularize the matrix before carrying out the computation. This can be faster for large matrices.

Example

julia> A0=[1 3; 4 5]; A1=[3 4; 5 6];
julia> id_op=S -> one(S) # Note: We use one(S) to be valid both for matrices and scalars
julia> exp_op=S -> exp(S)
julia> nep=SPMF_NEP([A0,A1],[id_op,exp_op]);
julia> compute_Mder(nep,1)-(A0+A1*exp(1))
2×2 Array{Float64,2}:
 0.0  0.0
 0.0  0.0

abstract  AbstractSPMF <: ProjectableNEP

An AbstractSPMF is an abstract class representing NEPs which can be represented as a sum of products of matrices and functions $M(λ)=Σ_i A_i f_i(λ)$, where i = 0,1,2,..., all of the matrices are of size $n×n$ and $f_i$ are functions.

Any AbstractSPMF has to have implementations of get_Av() and get_fv() which return the functions and matrices.

get_Av(nep::AbstractSPMF)

Returns an array of matrices $A_i$ in the AbstractSPMF: $M(λ)=Σ_i A_i f_i(λ)$

get_Av(nep::AbstractSPMF)

Returns an Array of functions (that can be evaluated both as scalar and matrix functions) $f_i$ in the AbstractSPMF: $M(λ)=Σ_i A_i f_i(λ)$

SumNEP{nep1::NEP,nep2::NEP}
SumNEP{nep1::AbstractSPMF,nep2::AbstractSPMF}

SumNEP is a function creating an object that corresponds to a sum of two NEPs, i.e., if nep is created by SumNEP it is defined by

where $M_1$ and $M_2$ are defined by nep1 and nep2.

Example:

julia> nep1=DEP([ones(3,3),randn(3,3)])
julia> nep2=PEP([ones(3,3),randn(3,3),randn(3,3)])
julia> sumnep=SumNEP(nep1,nep2);
julia> s=3.0;
julia> M=compute_Mder(sumnep,s);
3×3 Array{Float64,2}:
  8.54014     6.71897   7.12007
 -0.943908  -13.0795   -0.621659
  6.03155    -7.26726  -6.42828
julia> M1=compute_Mder(nep1,s);
julia> M2=compute_Mder(nep2,s);
julia> M1+M2  # Same as M
3×3 Array{Float64,2}:
  8.54014     6.71897   7.12007
 -0.943908  -13.0795   -0.621659
  6.03155    -7.26726  -6.42828

See also: SPMFSumNEP, GenericSumNEP

struct GenericSumNEP{NEP1<:NEP,NEP2<:NEP}  <: NEP

See also: SumNEP, SPMFSumNEP

struct SPMFSumNEP{NEP1<:AbstractSPMF,NEP2<:AbstractSPMF}  <: AbstractSPMF{AbstractMatrix}

See also: SumNEP, GenericSumNEP

struct LowRankFactorizedNEP <: AbstractSPMF
function LowRankFactorizedNEP(L::Vector,U::Vector,f::Vector)
function LowRankFactorizedNEP(L::Vector,U::Vector,A::Vector, f::Vector)

Representation of a NEP which has low rank in the sense that it is an SPMF where each of the terms are factorized: A[i]=L[i]*U[i]'. The factorization is provided in the L and U vectors and the full matrix A[i] can be either provided (or is otherwise implicitly computed).

Example:

julia> L=randn(5,1); U=randn(5,1);
julia> f=S->exp(S)
julia> nep=LowRankFactorizedNEP([L],[U],[f]);
julia> X=randn(5,2);
julia> norm(compute_Mlincomb(nep,0.0,X)-L*U'*X*ones(2),1)
6.661338147750939e-16

struct CORKPencil
function CORKPencil(M,N,Av,Bv,Z);
function CORKPencil(nep,is)

The struct CORKPencil represents a pencil with a particular structure, as given in the reference. The data can either be constructed directly via the first constructor, or from a NEP in the second constructor. The second constructor takes a NEP and is which specifies a CORK-structure as well as an approximation method. This can be objects of the type IarCorkLinearization or NleigsCorkLinearization, which are the CORK-linearizations equivalent to (certain versions of) iar and nleigs.

See buildPencil how to build standard pencil.

Example:

The following example constructs a CORKPencil from a general NEP and then computes approximations of NEPs by the interpolation approach of nleigs.

julia> using LinearAlgebra
julia> A=(1:4)*(1:4)'+I; B=diagm(1 => [1,2,3]); C=ones(4,4);
julia> f1= λ-> one(λ);
julia> f2= λ-> λ;
julia> f3= λ-> exp(sin(λ/2));
julia> nep=SPMF_NEP([A,B,C],[f1,f2,f3]);
julia> cp=CORKPencil(nep,NleigsCorkLinearization());
julia> (A,B)=buildPencil(cp) # Form the pencil
julia> λv=eigen(A,B).values;
julia> λ=λv[sortperm(abs.(λv))[1]]; # Smallest eigval
julia> minimum(svdvals(compute_Mder(nep,λ))) # It's a solution
2.4364382475487156e-11

References:

• Compact rational Krylov methods for nonlinear eigenvalue problems SIAM Journal on Matrix Analysis and Applications, 36 (2), 820-838, 2015.

(A,B)=buildPencil(cp)

Constructs a pencil from a CORKPencil or CORKPencilLR. The returned matrices correspond to a generalized eigenvalue problem. See also CORKPencil, CORKPencilLR. See [lowRankCompress] for examples.

struct CORKPencilLR
function CORKPencilLR(M,N,Av,AvLR,Bv,BvLR,Z);

Represents / constructs a low-rank CORK-pencil. AvLR, BvLR and Z correspond to the low-rank factorization of terms Av and Bv. See CORKPencil and reference below.

Example

The example illustrate a low-rank linearization of a Taylor expansion of the NEP $A0-λI+vv^Te^{-λ}$.

julia> A0=[1.0 3.0; -1.0 2.0]/10;
julia> v=reshape([-1.0 ; 1]/sqrt(2),n,1);

julia> Av=[-A0-v*v']
julia> Bv=[-one(A0)-v*v']
julia> BvLR=[v/2, -v/3, v/4, -v/5, v/6, -v/7,  v/8, -v/9]
julia> AvLR=zero.(BvLR);
julia> Z=v;
julia> d=9;
julia> M=diagm( 0 =>  ones(d) )[2:end,:]
julia> N=diagm( -1 =>  1 ./ (1:d-1) )[2:end,:]
julia> cplr=CORKPencilLR(M,N,Av,AvLR,Bv,BvLR,Z);
julia> (AA,BB)=buildPencil(cplr);
julia> λ=eigen(AA,BB).values[end];
julia> minimum(svdvals(A0-λ*I+v*v'*exp(-λ)))
8.870165379112754e-13

References:

• Section 7 in Compact rational Krylov methods for nonlinear eigenvalue problems SIAM Journal on Matrix Analysis and Applications, 36 (2), 820-838, 2015.

cplr=lowRankCompress(cp_org::CORKPencil,dtilde,rk)

Constructs a CORKPencilLR from a CORKPencil. This is done by assuming that terms higher than dtilde are of low rank, with rank rk. More precisely, all A[j] and B[j] for j>dtilde are assumed to be of the form $C_jZ^T$.

Example

This illustrates how to form a CORKPencil from a NEP, and subsequently form a smaller pencil using CORKPencilLR.

julia> A0=[1.0 3.0; -1.0 2.0]/10;
julia> v=[-1.0 ; 1]/sqrt(2);
julia> nep=DEP([A0,v*v']);
julia> cp_org=CORKPencil(nep,IarCorkLinearization(d=10));
julia> cplr=lowRankCompress(cp_org::CORKPencil,1,1);
julia> (AA,BB)=buildPencil(cplr);
julia> λ=eigen(AA,BB).values[end];
julia> minimum(svdvals(A0-λ*I+v*v'*exp(-λ))) # Check if it is an eigval
2.7621446071952216e-14