KrylovKit.jl

A KrylovFactorization represents a partial factorization or decomposition of a linear map.

Concrete Implementations:

• LanczosFactorization: For symmetric/Hermitian matrices.
• BlockLanczosFactorization: For block versions of Lanczos.
• ArnoldiFactorization: For general matrices.
• GKLFactorization: For Golub-Kahan-Lanczos bidiagonalization.

Accessing Components: You can access the defining components of a KrylovFactorization or GKLFactorization using these functions:

• basis(f): Returns the basis.
• rayleighquotient(f): Returns the Rayleigh quotient.
• residual(f): Returns the residual.
• normres(f): Returns the norm of the residual.
• rayleighextension(f): Returns the Rayleigh extension (often a SimpleBasisVector).

Specialized Types:

• SimpleBasisVector: A dedicated type to represent the Rayleigh extension without allocating a full vector.
• PackedHessenberg: A custom matrix type used to store the Hessenberg matrix of an Arnoldi factorization in a packed format (omitting zeros).

KrylovKit uses the Basis{T} abstract type to represent a basis of vectors (e.g., a Krylov subspace). For algorithms requiring an orthogonal basis, use the OrthonormalBasis{T} implementation.

To manipulate these bases, you can use the following operations:

• Orthogonalization: Use orthogonalize or orthonormalize to orthogonalize/orthonormalize a vector against another vector or an existing OrthonormalBasis. Use the in-place versions orthogonalize!! and orthonormalize!! to modify vectors in place.
• Projection: Use project!! to obtain the expansion coefficients of a vector in terms of an orthonormal basis. Use unproject!! to perform the inverse calculation.
• Basis Updates: Use rank1update! to transform an orthonormal basis via a rank-1 update (this changes the subspace). Use basistransform! to perform a rotation of the basis that does not change the spanned subspace.

# Example conceptual usage
basis = OrthonormalBasis(T, vectors)
# Orthogonalize a new vector against the basis
orthonormalize!(basis, new_vector)

# Project a vector onto the basis to get coefficients
coeffs = project!!(basis, vector)

High-level functions such as linsolve, eigsolve, geneigsolve, svdsolve, exponentiate, and expintegrator follow a consistent interface.

Input: The Linear Map

The first argument A represents the linear map. This can be:

1. An instance of AbstractMatrix.
2. A function or callable object that encodes the action of the linear map on a vector.
3. A Julia do block that implements the action on a vector x.

Input: Arguments and Keywords

• args...: Additional arguments required by the specific solver.
• kwargs...: Configuration keywords:
  • Linear map properties: issymmetric, ishermitian, isposdef.
  • Solution strategy: tol (tolerance), krylovdim (Krylov subspace dimension), maxiter (maximum iterations).
  • Output control: verbosity.

Output: Results and Convergence Info

Functions return one or more solution entries followed by a final info object of type ConvergeInfo. This object contains convergence status, residuals, the norm of the residual, and the number of operations used.

Verbosity Levels

Control how much information is printed to STDOUT using the verbosity keyword:

• SILENT_LEVEL (default): No information printed.
• WARN_LEVEL: A single message at the end (warning if failed, info if succeeded).
• STARTSTOP_LEVEL: Information about the current state is displayed after every iteration.
• verbosity > STARTSTOP_LEVEL: Detailed information about individual Krylov expansion steps is displayed.

# Using a do-block to define the linear map
results..., info = linsolve(args...; kwargs...) do x
    y = # implement linear map on x
    return y
end

The eigsolve and realeigsolve routines are compatible with reverse-mode automatic differentiation via the ChainRules ecosystem.

Key Details:

• Adjoint Map: The adjoint problem requires the adjoint of the linear map. If the map is an AbstractMatrix, its adjoint is used. If it is a function f, the AD engine computes the adjoint via ChainRulesCore.rrule_via_ad.
• Selecting Adjoint Algorithms: You can select how the adjoint problem is solved using the alg_rrule keyword argument:
  • If a linear solver (e.g., GMRES, BiCGStab) is specified, the adjoint problem solves multiple linear problems.
  • If an eigenvalue solver (e.g., Arnoldi) is specified, the adjoint problem is solved as a single, larger eigenvalue problem.
• Gauge Invariance Warning: Because eigenvector phase is not unique, cost functions depending on eigenvectors should be 'gauge invariant'. If a cost function is 'gauge dependent', KrylovKit will print a warning. To suppress this warning, pass an algorithm to alg_rrule with verbosity=SILENT_LEVEL-1.

KrylovKit.jl is designed for flexibility regarding the types of inputs it accepts:

Linear Maps

KrylovKit accepts general functions or callable objects as linear maps.

• If the first argument (the linear map) is a subtype of AbstractMatrix, the package uses matrix-vector multiplication.
• Otherwise, the map is applied as a function call.
• You do not need to wrap your function in LinearMap or LinearOperator types.

Vectors

KrylovKit does not require vectors to be subtypes of AbstractVector. It supports any Julia object that follows the VectorInterface.jl interface. This includes:

• Higher-dimensional arrays.
• Custom user types.
• Tuples and NamedTuples.
• Arbitrarily nested combinations of tuples and arrays.
• CuArray objects (for GPU execution of the linear map application).

Note: While the linear map application can run on a GPU (e.g., using CuArray), the internal Krylov subspace computations (like diagonalizing the projected matrix) are always performed on the CPU.

The svdsolve routine is compatible with reverse-mode automatic differentiation (AD) via the ChainRules ecosystem.

To use AD, the linear map $A$ must support both its action and its adjoint (conjugate transpose). The adjoint problem can be solved using different implementations via the alg_rrule keyword argument:

1. Linear Solver Approach: If you specify a linear solver like GMRES or BiCGStab for alg_rrule, the adjoint problem is solved as a number of linear problems equal to the number of requested singular values/vectors.
2. Eigenvalue Solver Approach: If an eigenvalue solver (such as Arnoldi) is specified, the adjoint problem is solved as a single, larger eigenvalue problem.

Important: Gauge Invariance Warning Singular vectors have an arbitrary phase freedom. If your cost function is 'gauge dependent' (meaning its value changes when left and right singular vectors are simultaneously changed by a common phase factor), KrylovKit will print a warning during the adjoint calculation.

To suppress this warning, pass an algorithm to alg_rrule with verbosity=SILENT_LEVEL-1.

# Example concept: selecting an rrule implementation
svdsolve(A, k, alg_rrule=GMRES(verbosity=SILENT_LEVEL-1))