Details of the implementation

Details of the implementation

Orthogonalization

To denote a basis of vectors, e.g. to represent a given Krylov subspace, there is an abstract type Basis{T}

KrylovKit.BasisType.
abstract type Basis{T} end

An abstract type to collect specific types for representing a basis of vectors of type T.

Implementations of Basis{T} behave in many ways like Vector{T} and should have a length, can be indexed (getindex and setindex!), iterated over (iterate), and support resizing (resize!, pop!, push!, empty!, sizehint!).

The type T denotes the type of the elements stored in an Basis{T} and can be any custom type that has vector like behavior (as defined in the docs of KrylovKit).

See OrthonormalBasis.

source

Many Krylov based algorithms use an orthogonal basis to parametrize the Krylov subspace. In that case, the specific implementation OrthonormalBasis{T} can be used:

OrthonormalBasis{T} <: Basis{T}

A list of vector like objects of type T that are mutually orthogonal and normalized to one, representing an orthonormal basis for some subspace (typically a Krylov subspace). See also Basis

Orthonormality of the vectors contained in an instance b of OrthonormalBasis (i.e. all(dot(b[i],b[j]) == I[i,j] for i=1:lenght(b), j=1:length(b))) is not checked when elements are added; it is up to the algorithm that constructs b to guarantee orthonormality.

One can easily orthogonalize or orthonormalize a given vector v with respect to a b::OrthonormalBasis using the functions w, = orthogonalize(v,b,...) or w, = orthonormalize(v,b,...). The resulting vector w of the latter can then be added to b using push!(b, w). Note that in place versions orthogonalize!(v, b, ...) or orthonormalize!(v, b, ...) are also available.

Finally, a linear combination of the vectors in b::OrthonormalBasis can be obtained by multiplying b with a Vector{<:Number} using * or mul! (if the output vector is already allocated).

source

We can orthogonalize or orthonormalize a given vector to another vector (assumed normalized) or to a given OrthonormalBasis.

orthogonalize(v, b::OrthonormalBasis, [x::AbstractVector,] algorithm::Orthogonalizer]) -> w, x
orthogonalize!(v, b::OrthonormalBasis, [x::AbstractVector,] algorithm::Orthogonalizer]) -> w, x

orthogonalize(v, q, algorithm::Orthogonalizer]) -> w, s
orthogonalize!(v, q, algorithm::Orthogonalizer]) -> w, s

Orthogonalize vector v against all the vectors in the orthonormal basis b using the orthogonalization algorithm algorithm, and return the resulting vector w and the overlap coefficients x of v with the basis vectors in b.

In case of orthogonalize!, the vector v is mutated in place. In both functions, storage for the overlap coefficients x can be provided as optional argument x::AbstractVector with length(x) >= length(b).

One can also orthogonalize v against a given vector q (assumed to be normalized), in which case the orthogonal vector w and the inner product s between v and q are returned.

Note that w is not normalized, see also orthonormalize.

For algorithms, see ClassicalGramSchmidt, ModifiedGramSchmidt, ClassicalGramSchmidt2, ModifiedGramSchmidt2, ClassicalGramSchmidtIR and ModifiedGramSchmidtIR.

source
orthonormalize(v, b::OrthonormalBasis, [x::AbstractVector,] algorithm::Orthogonalizer]) -> w, β, x
orthonormalize!(v, b::OrthonormalBasis, [x::AbstractVector,] algorithm::Orthogonalizer]) -> w, β, x

orthonormalize(v, q, algorithm::Orthogonalizer]) -> w, β, s
orthonormalize!(v, q, algorithm::Orthogonalizer]) -> w, β, s

Orthonormalize vector v against all the vectors in the orthonormal basis b using the orthogonalization algorithm algorithm, and return the resulting vector w (of norm 1), its norm β after orthogonalizing and the overlap coefficients x of v with the basis vectors in b, such that v = β * w + b * x.

In case of orthogonalize!, the vector v is mutated in place. In both functions, storage for the overlap coefficients x can be provided as optional argument x::AbstractVector with length(x) >= length(b).

One can also orthonormalize v against a given vector q (assumed to be normalized), in which case the orthonormal vector w, its norm β before normalizing and the inner product s between v and q are returned.

See orthogonalize if w does not need to be normalized.

For algorithms, see ClassicalGramSchmidt, ModifiedGramSchmidt, ClassicalGramSchmidt2, ModifiedGramSchmidt2, ClassicalGramSchmidtIR and ModifiedGramSchmidtIR.

source

Dense linear algebra

KrylovKit relies on Julia's LinearAlgebra module from the standard library for most of its dense linear algebra dependencies.

Krylov factorizations

The central ingredient in a Krylov based algorithm is a Krylov factorization or decomposition of a linear map. Such partial factorizations are represented as a KrylovFactorization, of which LanczosFactorization and ArnoldiFactorization are two concrete implementations:

abstract type KrylovFactorization{T,S<:Number}
mutable struct LanczosFactorization{T,S<:Real}    <: KrylovFactorization{T,S}
mutable struct ArnoldiFactorization{T,S<:Number}  <: KrylovFactorization{T,S}

Structures to store a Krylov factorization of a linear map A of the form

    A * V = V * B + r * b'.

For a given Krylov factorization fact of length k = length(fact), the basis A is obtained via [basis(fact)](@ref basis) and is an instance of some subtype of [Basis{T}](@ref Basis), with alsolength(V) == kand whereTdenotes the type of vector like objects used in the problem. The Rayleigh quotientBis obtained as [rayleighquotient(fact)](@ref) andtypeof(B)is some subtype ofAbstractMatrix{S}withsize(B) == (k,k), typically a structured matrix. The residualris obtained as [residual(fact)](@ref) and is of typeT. One can also query [normres(fact)](@ref) to obtainnorm(r), the norm of the residual. The vectorbhas no dedicated name and often takes a default form (see below). It should be a subtype ofAbstractVectorof lengthkand can be obtained as [rayleighextension(fact)`](@ref) (by lack of a better dedicated name).

In particular, LanczosFactorization stores a Lanczos factorization of a real symmetric or complex hermitian linear map and has V::OrthonormalBasis{T} and B::SymTridiagonal{S<:Real}. ArnoldiFactorization stores an Arnoldi factorization of a general linear map and has V::OrthonormalBasis{T} and B::PackedHessenberg{S<:Number}. In both cases, b takes the default value $e_k$, i.e. the unit vector of all zeros and a one in the last entry, which is represented using SimpleBasisVector.

A Krylov factorization fact can be destructured as V, B, r, nr, b = fact with nr = norm(r).

LanczosFactorization and ArnoldiFactorization are mutable because they can expand! or shrink!. See also KrylovIterator (and in particular LanczosIterator and ArnoldiIterator) for iterators that construct progressively expanding Krylov factorizations of a given linear map and a starting vector.

source

A KrylovFactorization can be destructered into its defining components using iteration, but these can also be accessed using the following functions

KrylovKit.basisFunction.
    basis(fact::KrylovFactorization)

Return the list of basis vectors of a KrylovFactorization, which span the Krylov subspace. The return type is a subtype of Basis{T}, where T represents the type of the vectors used by the problem.

source
rayleighquotient(fact::KrylovFactorization)

Return the Rayleigh quotient of a KrylovFactorization, i.e. the reduced matrix within the basis of the Krylov subspace. The return type is a subtype of AbstractMatrix{<:Number}, typically some structured matrix type.

source
KrylovKit.residualFunction.
residual(fact::KrylovFactorization)

Return the residual of a KrylovFactorization. The return type is some vector of the same type as used in the problem. See also normres(F) for its norm, which typically has been computed already.

source
KrylovKit.normresFunction.
normres(fact::KrylovFactorization)

Return the norm of the residual of a KrylovFactorization. As this has typically already been computed, it is cheaper than (but otherwise equivalent to) norm(residual(F)).

source
rayleighextension(fact::KrylovFactorization)

Return the vector b appearing in the definition of a KrylovFactorization.

source

Krylov iterators

Given a linear map $A$ and a starting vector $x₀$, a Krylov factorization is obtained by sequentially building a Krylov subspace ${x₀, A x₀, A² x₀, ...}$. Rather then using this set of vectors as a basis, an orthonormal basis is generated by a process known as Lanczos or Arnoldi iteration (for symmetric/hermitian and for general matrices, respectively). These processes are represented as iterators in Julia:

abstract type KrylovIterator{F,T}
struct LanczosIterator{F,T,O<:Orthogonalizer} <: KrylovIterator{F,T}
struct ArnoldiIterator{F,T,O<:Orthogonalizer} <: KrylovIterator{F,T}

LanczosIterator(f, v₀, [orth::Orthogonalizer = KrylovDefaults.orth], keepvecs::Bool = true)
ArnoldiIterator(f, v₀, [orth::Orthogonalizer = KrylovDefaults.orth])

Iterators that take a linear map of type F and an initial vector of type T and generate an expanding KrylovFactorization thereof.

In particular, for a real symmetric or complex hermitian linear map f, LanczosIterator uses the Lanczos iteration scheme to build a successively expanding LanczosFactorization. While f cannot be tested to be symmetric or hermitian directly when the linear map is encoded as a general callable object or function, it is tested whether the imaginary part of dot(v, f(v)) is sufficiently small to be neglected.

Similarly, for a general linear map f, ArnoldiIterator iterates over progressively expanding ArnoldiFactorizations using the Arnoldi iteration.

The optional argument orth specifies which Orthogonalizer to be used. The default value in KrylovDefaults is to use ModifiedGramSchmidtIR, which possibly uses reorthogonalization steps. For LanczosIterator, one can use to discard the old vectors that span the Krylov subspace by setting the final argument keepvecs to false. This, however, is only possible if an orth algorithm is used that does not rely on reorthogonalization, such as ClassicalGramSchmidt() or ModifiedGramSchmidt(). In that case, the iterator strictly uses the Lanczos three-term recurrence relation.

When iterating over an instance of KrylovIterator, the values being generated are subtypes of KrylovFactorization, which can be immediately destructured into a basis, rayleighquotient, residual, normres and rayleighextension, for example as

for V,B,r,nr,b in ArnoldiIterator(f, v₀)
    # do something
    nr < tol && break # a typical stopping criterion
end

Note, however, that if keepvecs=false in LanczosIterator, the basis V cannot be extracted. Since the iterators don't know the dimension of the underlying vector space of objects of type T, they keep expanding the Krylov subspace until normres falls below machine precision eps for the given eltype(T).

The internal state of LanczosIterator and ArnoldiIterator is the same as the return value, i.e. the corresponding LanczosFactorization or ArnoldiFactorization. However, as Julia's Base iteration interface (using Base.iterate) requires that the state is not mutated, a deepcopy is produced upon every next iteration step.

Instead, you can also mutate the KrylovFactorization in place, using the following interface, e.g. for the same example above

iterator = ArnoldiIterator(f, v₀)
factorization = initialize(iterator)
while normres(factorization) > tol
    expand!(iterator, f)
    V,B,r,nr,b = factorization
    # do something
end

Here, initialize(::KrylovIterator) produces the first Krylov factorization of length 1, and expand!(::KrylovIterator,::KrylovFactorization)(@ref) expands the factorization in place. See also initialize!(::KrylovIterator,::KrylovFactorization) to initialize in an already existing factorization (most information will be discarded) and shrink!(::KrylovFactorization, k) to shrink an existing factorization down to length k.

source

The following functions allow to manipulate a KrylovFactorization obtained from such a KrylovIterator:

KrylovKit.expand!Function.
expand!(iter::KrylovIteraotr, fact::KrylovFactorization)

Expand the Krylov factorization fact by one using the linear map and parameters in iter.

source
KrylovKit.shrink!Function.
shrink!(fact::KrylovFactorization, k)

Shrink an existing Krylov factorization fact down to have length k. Does nothing if length(fact)<=k.

source
KrylovKit.initialize!Function.
initialize!(iter::KrylovIteraotr, fact::KrylovFactorization)

Initialize a length 1 Kryov factorization corresponding to iter in the already existing factorization fact, thereby destroying all the information it currently holds.

source
KrylovKit.initializeFunction.
initialize(iter::KrylovIteraotr)

Initialize a length 1 Kryov factorization corresponding to iter.

source