# Implementation details

## 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 for a specific implementation.

source

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

KrylovKit.OrthonormalBasisType
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:length(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 KrylovKit.OrthonormalBasis.

KrylovKit.orthogonalizeFunction
orthogonalize(v, b::OrthonormalBasis, [x::AbstractVector,] alg::Orthogonalizer]) -> w, x
orthogonalize!(v, b::OrthonormalBasis, [x::AbstractVector,] alg::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 alg of type Orthogonalizer, 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 more information on possible orthogonalization algorithms, see Orthogonalizer and its concrete subtypes ClassicalGramSchmidt, ModifiedGramSchmidt, ClassicalGramSchmidt2, ModifiedGramSchmidt2, ClassicalGramSchmidtIR and ModifiedGramSchmidtIR.

source
KrylovKit.orthonormalizeFunction
orthonormalize(v, b::OrthonormalBasis, [x::AbstractVector,] alg::Orthogonalizer]) -> w, β, x
orthonormalize!(v, b::OrthonormalBasis, [x::AbstractVector,] alg::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 alg of type Orthogonalizer, 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 more information on possible orthogonalization algorithms, see Orthogonalizer and its concrete subtypes ClassicalGramSchmidt, ModifiedGramSchmidt, ClassicalGramSchmidt2, ModifiedGramSchmidt2, ClassicalGramSchmidtIR and ModifiedGramSchmidtIR.

source

The expansion coefficients of a general vector in terms of a given orthonormal basis can be obtained as

KrylovKit.project!Function
project!(y::AbstractVector, b::OrthonormalBasis, x,
[α::Number = 1, β::Number = 0, r = Base.OneTo(length(b))])

For a given orthonormal basis b, compute the expansion coefficients y resulting from projecting the vector x onto the subspace spanned by b; more specifically this computes

    y[j] = β*y[j] + α * dot(b[r[j]], x)

for all $j ∈ r$.

source

whereas the inverse calculation is obtained as

KrylovKit.unproject!Function
unproject!(y, b::OrthonormalBasis, x::AbstractVector,
[α::Number = 1, β::Number = 0, r = Base.OneTo(length(b))])

For a given orthonormal basis b, reconstruct the vector-like object y that is defined by expansion coefficients with respect to the basis vectors in b in x; more specifically this computes

    y = β*y + α * sum(b[r[i]]*x[i] for i = 1:length(r))
source

An orthonormal basis can be transformed using a rank-1 update using

KrylovKit.rank1update!Function
rank1update!(b::OrthonormalBasis, y, x::AbstractVector,
[α::Number = 1, β::Number = 1, r = Base.OneTo(length(b))])

Perform a rank 1 update of a basis b, i.e. update the basis vectors as

    b[r[i]] = β*b[r[i]] + α * y * conj(x[i])

It is the user's responsibility to make sure that the result is still an orthonormal basis.

source

Note that this changes the subspace. A mere rotation of the basis, which does not change the subspace spanned by it, can be computed using

KrylovKit.basistransform!Function
basistransform!(b::OrthonormalBasis, U::AbstractMatrix)

Transform the orthonormal basis b by the matrix U. For b an orthonormal basis, the matrix U should be real orthogonal or complex unitary; it is up to the user to ensure this condition is satisfied. The new basis vectors are given by

    b[j] ← b[i] * U[i,j]

and are stored in b, so the old basis vectors are thrown away. Note that, by definition, the subspace spanned by these basis vectors is exactly the same.

source

## Dense linear algebra

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

## Factorization types

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.KrylovFactorizationType
abstract type KrylovFactorization{T,S<:Number}

Abstract type 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 V is obtained via basis(fact) and is an instance of some subtype of Basis{T}, with also length(V) == k and where T denotes the type of vector like objects used in the problem. The Rayleigh quotient B is obtained as rayleighquotient(fact) and typeof(B) is some subtype of AbstractMatrix{S} with size(B) == (k,k), typically a structured matrix. The residual r is obtained as residual(fact) and is of type T. One can also query normres(fact) to obtain norm(r), the norm of the residual. The vector b has no dedicated name and often takes a default form (see below). It should be a subtype of AbstractVector of length k and can be obtained as rayleighextension(fact) (by lack of a better dedicated name).

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

See also LanczosFactorization and ArnoldiFactorization for concrete implementations, and KrylovIterator (with in particular LanczosIterator and ArnoldiIterator) for iterators that construct progressively expanding Krylov factorizations of a given linear map and a starting vector.

source
KrylovKit.LanczosFactorizationType
mutable struct LanczosFactorization{T,S<:Real} <: KrylovFactorization{T,S}

Structure to store a Lanczos factorization of a real symmetric or complex hermitian linear map A of the form

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

For a given Lanczos factorization fact of length k = length(fact), the basis V is obtained via basis(fact) and is an instance of OrthonormalBasis{T}, with also length(V) == k and where T denotes the type of vector like objects used in the problem. The Rayleigh quotient B is obtained as rayleighquotient(fact) and is of type SymTridiagonal{S<:Real} with size(B) == (k,k). The residual r is obtained as residual(fact) and is of type T. One can also query normres(fact) to obtain norm(r), the norm of the residual. The vector b has no dedicated name but can be obtained via rayleighextension(fact). It 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 Lanczos factorization fact can be destructured as V, B, r, nr, b = fact with nr = norm(r).

LanczosFactorization is mutable because it can expand! or shrink!. See also LanczosIterator for an iterator that constructs a progressively expanding Lanczos factorizations of a given linear map and a starting vector. See ArnoldiFactorization and ArnoldiIterator for a Krylov factorization that works for general (non-symmetric) linear maps.

source
KrylovKit.ArnoldiFactorizationType
mutable struct ArnoldiFactorization{T,S} <: KrylovFactorization{T,S}

Structure to store an Arnoldi factorization of a linear map A of the form

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

For a given Arnoldi factorization fact of length k = length(fact), the basis V is obtained via basis(fact) and is an instance of OrthonormalBasis{T}, with also length(V) == k and where T denotes the type of vector like objects used in the problem. The Rayleigh quotient B is obtained as rayleighquotient(fact) and is of type B::PackedHessenberg{S<:Number} with size(B) == (k,k). The residual r is obtained as residual(fact) and is of type T. One can also query normres(fact) to obtain norm(r), the norm of the residual. The vector b has no dedicated name but can be obtained via rayleighextension(fact). It 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.

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

ArnoldiFactorization is mutable because it can expand! or shrink!. See also ArnoldiIterator for an iterator that constructs a progressively expanding Arnoldi factorizations of a given linear map and a starting vector. See LanczosFactorization and LanczosIterator for a Krylov factorization that is optimized for real symmetric or complex hermitian linear maps.

source
KrylovKit.GKLFactorizationType
mutable struct GKLFactorization{TU,TV,S<:Real}

Structure to store a Golub-Kahan-Lanczos (GKL) bidiagonal factorization of a linear map A of the form

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

For a given GKL factorization fact of length k = length(fact), the two bases U and V are obtained via basis(fact, :U) and basis(fact, :V). Here, U and V are instances of OrthonormalBasis{T}, with also length(U) == length(V) == k and where T denotes the type of vector like objects used in the problem. The Rayleigh quotient B is obtained as rayleighquotient(fact) and is of type Bidiagonal{S<:Number} with size(B) == (k,k). The residual r is obtained as residual(fact) and is of type T. One can also query normres(fact) to obtain norm(r), the norm of the residual. The vector b has no dedicated name but can be obtained via rayleighextension(fact). It 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 GKL factorization fact can be destructured as U, V, B, r, nr, b = fact with nr = norm(r).

GKLFactorization is mutable because it can expand! or shrink!. See also GKLIterator for an iterator that constructs a progressively expanding GKL factorizations of a given linear map and a starting vector u₀.

source

A KrylovFactorization or GKLFactorization can be destructured 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

As the rayleighextension is typically a simple basis vector, we have created a dedicated type to represent this without having to allocate an actual vector, i.e.

KrylovKit.SimpleBasisVectorType
SimpleBasisVector(m, k)

Construct a simple struct SimpleBasisVector <: AbstractVector{Bool} representing a coordinate basis vector of length m in the direction of k, i.e. for e_k = SimpleBasisVector(m, k) we have length(e_k) = m and e_k[i] = (i == k).

source

Furthermore, to store the Rayleigh quotient of the Arnoldi factorization in a manner that can easily be expanded, we have constructed a custom matrix type to store the Hessenberg matrix in a packed format (without zeros):

KrylovKit.PackedHessenbergType
struct PackedHessenberg{T,V<:AbstractVector{T}} <: AbstractMatrix{T}
data::V
n::Int
end

A custom struct to store a Hessenberg matrix in a packed format (without zeros). Hereto, the non-zero entries are stored sequentially in vector data of length n(n+1)/2.

source

## Factorization 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.KrylovIteratorType
abstract type KrylovIterator{F,T}

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

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.

See LanczosIterator and ArnoldiIterator for concrete implementations and more information.

source
KrylovKit.LanczosIteratorType
struct LanczosIterator{F,T,O<:Orthogonalizer} <: KrylovIterator{F,T}
LanczosIterator(f, v₀, [orth::Orthogonalizer = KrylovDefaults.orth, keepvecs::Bool = true])

Iterator that takes a linear map f::F (supposed to be real symmetric or complex hermitian) and an initial vector v₀::T and generates an expanding LanczosFactorization thereof. In particular, LanczosIterator uses the Lanczos iteration scheme to build a successively expanding Lanczos factorization. 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.

The argument f can be a matrix, or a function accepting a single argument v, so that f(v) implements the action of the linear map on the vector v.

The optional argument orth specifies which Orthogonalizer to be used. The default value in KrylovDefaults is to use ModifiedGramSchmidtIR, which possibly uses reorthogonalization steps. 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 LanczosIterator, the values being generated are instances of LanczosFactorization, which can be immediately destructured into a basis, rayleighquotient, residual, normres and rayleighextension, for example as

for (V, B, r, nr, b) in LanczosIterator(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 iterator does not know the dimension of the underlying vector space of objects of type T, it keeps expanding the Krylov subspace until the residual norm nr falls below machine precision eps(typeof(nr)).

The internal state of LanczosIterator is the same as the return value, i.e. the corresponding LanczosFactorization. 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 = LanczosIterator(f, v₀)
factorization = initialize(iterator)
while normres(factorization) > tol
expand!(iterator, factorization)
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
KrylovKit.ArnoldiIteratorType
struct ArnoldiIterator{F,T,O<:Orthogonalizer} <: KrylovIterator{F,T}
ArnoldiIterator(f, v₀, [orth::Orthogonalizer = KrylovDefaults.orth])

Iterator that takes a general linear map f::F and an initial vector v₀::T and generates an expanding ArnoldiFactorization thereof. In particular, ArnoldiIterator iterates over progressively expanding Arnoldi factorizations using the Arnoldi iteration.

The argument f can be a matrix, or a function accepting a single argument v, so that f(v) implements the action of the linear map on the vector v.

The optional argument orth specifies which Orthogonalizer to be used. The default value in KrylovDefaults is to use ModifiedGramSchmidtIR, which possibly uses reorthogonalization steps.

When iterating over an instance of ArnoldiIterator, the values being generated are instances of ArnoldiFactorization, 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

Since the iterator does not know the dimension of the underlying vector space of objects of type T, it keeps expanding the Krylov subspace until the residual norm nr falls below machine precision eps(typeof(nr)).

The internal state of ArnoldiIterator is the same as the return value, i.e. the corresponding 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 ArnoldiFactorization 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, factorization)
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

Similarly, there is also an iterator for the Golub-Kahan-Lanczos bidiagonalization proces:

KrylovKit.GKLIteratorType
struct GKLIterator{F,TU,O<:Orthogonalizer}
GKLIterator(f, u₀, [orth::Orthogonalizer = KrylovDefaults.orth, keepvecs::Bool = true])

Iterator that takes a general linear map f::F and an initial vector u₀::TU and generates an expanding GKLFactorization thereof. In particular, GKLIterator implements the Golub-Kahan-Lanczos bidiagonalization procedure. Note, however, that this implementation starts from a vector u₀ in the codomain of the linear map f, which will end up (after normalisation) as the first column of U.

The argument f can be a matrix, a tuple of two functions where the first represents the normal action and the second the adjoint action, or a function accepting two arguments, where the first argument is the vector to which the linear map needs to be applied, and the second argument is either Val(false) for the normal action and Val(true) for the adjoint action. Note that the flag is thus a Val type to allow for type stability in cases where the vectors in the domain and the codomain of the linear map have a different type.

The optional argument orth specifies which Orthogonalizer to be used. The default value in KrylovDefaults is to use ModifiedGramSchmidtIR, which possibly uses reorthogonalization steps.

When iterating over an instance of GKLIterator, the values being generated are instances fact of GKLFactorization, which can be immediately destructured into a basis(fact, :U), basis(fact, :V), rayleighquotient, residual, normres and rayleighextension, for example as

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

Since the iterator does not know the dimension of the underlying vector space of objects of type T, it keeps expanding the Krylov subspace until the residual norm nr falls below machine precision eps(typeof(nr)).

The internal state of GKLIterator is the same as the return value, i.e. the corresponding GKLFactorization. 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 GKLFactorization in place, using the following interface, e.g. for the same example above

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

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

source

As an alternative to the standard iteration interface from Julia Base (using iterate), these iterative processes and the factorizations they produce can also be manipulated using the following functions:

KrylovKit.expand!Function
expand!(iter::KrylovIterator, 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::KrylovIterator, fact::KrylovFactorization)

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

source