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.Basis — Type.abstract type Basis{T} endAn 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.
Many Krylov based algorithms use an orthogonal basis to parametrize the Krylov subspace. In that case, the specific implementation OrthonormalBasis{T} can be used:
KrylovKit.OrthonormalBasis — Type.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).
We can orthogonalize or orthonormalize a given vector to another vector (assumed normalized) or to a given OrthonormalBasis.
KrylovKit.orthogonalize — Function.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, sOrthogonalize 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.
KrylovKit.orthonormalize — Function.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, β, sOrthonormalize 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.
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:
KrylovKit.KrylovFactorization — Type.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.
A KrylovFactorization can be destructered into its defining components using iteration, but these can also be accessed using the following functions
KrylovKit.basis — Function. 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.
KrylovKit.rayleighquotient — Function.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.
KrylovKit.residual — Function.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.
KrylovKit.normres — Function.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)).
KrylovKit.rayleighextension — Function.rayleighextension(fact::KrylovFactorization)Return the vector b appearing in the definition of a KrylovFactorization.
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:
KrylovKit.KrylovIterator — Type.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
endNote, 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
endHere, 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.
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.
KrylovKit.shrink! — Function.shrink!(fact::KrylovFactorization, k)Shrink an existing Krylov factorization fact down to have length k. Does nothing if length(fact)<=k.
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.
KrylovKit.initialize — Function.initialize(iter::KrylovIteraotr)Initialize a length 1 Kryov factorization corresponding to iter.