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} 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
.
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, 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
.
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, β, 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
.
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 also
length(V) == kand where
Tdenotes the type of vector like objects used in the problem. The Rayleigh quotient
Bis obtained as [
rayleighquotient(fact)](@ref) and
typeof(B)is some subtype of
AbstractMatrix{S}with
size(B) == (k,k), typically a structured matrix. The residual
ris obtained as [
residual(fact)](@ref) and is of type
T. One can also query [
normres(fact)](@ref) to obtain
norm(r), the norm of the residual. The vector
bhas no dedicated name and often takes a default form (see below). It should be a subtype of
AbstractVectorof length
kand 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
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
.
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
.