Implementation details

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.

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(inner(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).

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.

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.

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.

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.

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] + α * inner(b[r[j]], x)

for all $j ∈ r$.

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))

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.

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.

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.

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.

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.

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₀.

    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.

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.

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.

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)).

rayleighextension(fact::KrylovFactorization)

Return the vector b appearing in the definition of a KrylovFactorization; often it is simply the last coordinate unit vector, which can be represented using SimpleBasisVector.

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).

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.

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.

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 inner(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.

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.

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.

expand!(iter::KrylovIterator, fact::KrylovFactorization)

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

shrink!(fact::KrylovFactorization, k)

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

initialize(iter::KrylovIterator)

Initialize a length 1 Krylov factorization corresponding to iter.

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.