Eigsolve initialization

Project.toml

@@ -12,9 +12,11 @@ VectorInterface = "409d34a3-91d5-4945-b6ec-7529ddf182d8"
 
 [weakdeps]
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [extensions]
 KrylovKitChainRulesCoreExt = "ChainRulesCore"
+KrylovKitSparseArraysExt = "SparseArrays"
 
 [compat]
 Aqua = "0.6, 0.7, 0.8"
@@ -26,6 +28,7 @@ Logging = "1"
 PackageExtensionCompat = "1"
 Printf = "1"
 Random = "1"
+SparseArrays = "1"
 Test = "1"
 TestExtras = "0.2,0.3"
 VectorInterface = "0.5"
@@ -38,9 +41,11 @@ ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
 ChainRulesTestUtils = "cdddcdb0-9152-4a09-a978-84456f9df70a"
 FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TestExtras = "5ed8adda-3752-4e41-b88a-e8b09835ee3a"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [targets]
-test = ["Test", "Aqua", "Logging", "TestExtras", "ChainRulesTestUtils", "ChainRulesCore", "FiniteDifferences", "Zygote"]
+test = ["Test", "Aqua", "Logging", "TestExtras", "SparseArrays", "ChainRulesTestUtils",
+        "ChainRulesCore", "FiniteDifferences", "Zygote"]

docs/src/man/implementation.md

@@ -1,5 +1,43 @@
 # Implementation details
 
+## Linear map interface
+
+KrylovKit.jl aims to support a large variety of linear maps. By default, you can use any
+`AbstractMatrix`, in which case `A * x` is used to compute the action of the linear map on
+a vector `x`. Alternatively, you can use any callable object, for which the action is
+encoded as `A(x)`. For custom objects that do not which to implement the callable interface,
+you may finally implement `KrylovKit.apply(A, x)` to specify the action.
+
+```@docs
+KrylovKit.apply
+```
+
+For some algorithms, the adjoint of the linear map is also required. In these cases, again
+the `AbstractMatrix` inputs will use `A * x` and `A' * x`. For callable objects, the
+required signature is `f(x, ::Val{true})` for the adjoint, and `f(x, ::Val{false})` for the
+regular action. Alternatively, you may specify a tuple of callable objects `(f, fadjoint)`,
+which will be used for the action and adjoint, respectively. Finally, you may implement
+`KrylovKit.apply_normal(A, x)` and `KrylovKit.apply_adjoint(A, x)` for the action and adjoint.
+
+```@docs
+KrylovKit.apply_normal
+KrylovKit.apply_adjoint
+```
+
+## Initial vector interface
+
+Some algorithms require an initial vector to start the Krylov expansion. This is best
+specified directly as an argument `x0`, which often serves as an initial guess. For
+convenience, linear maps that are `AbstractMatrix` will generate an initial vector through
+`Random.rand!(similar(A, scalartype(A), size(A, 1)))`. For function handles it is usually
+not possible to deduce the type of the initial vector, so you have to specify the initial
+vector explicitly. For specific custom structs, you can enable support for generating
+initial vectors by extending [`initialize_vector`](@ref).
+
+```@docs
+KrylovKit.initialize_vector
+```
+
 ## Orthogonalization
 To denote a basis of vectors, e.g. to represent a given Krylov subspace, there is an
 abstract type `Basis{T}`

docs/src/man/intro.md

@@ -58,7 +58,8 @@ is a keyword argument `verbosity` that determines how much information is printe
 `verbosity = 1`, a single message at the end of the algorithm will be displayed, which is a
 warning if the algorithm did not succeed in finding the solution, or some information if it
 did. For `verbosity = 2`, information about the current state is displayed after every
-iteration of the algorithm. Finally, for `verbosity > 2`, information about the individual Krylov expansion steps is displayed.
+iteration of the algorithm. Finally, for `verbosity > 2`, information about the individual
+Krylov expansion steps is displayed.
 
 The return value contains one or more entries that define the solution, and a final
 entry `info` of type `ConvergeInfo` that encodes information about the solution, i.e.

ext/KrylovKitSparseArraysExt.jl

@@ -0,0 +1,9 @@
+module KrylovKitSparseArraysExt
+
+using SparseArrays
+using KrylovKit
+using VectorInterface: scalartype
+
+KrylovKit.initialize_vector(A::SparseMatrixCSC) = rand(scalartype(A), size(A, 1))
+
+end

src/KrylovKit.jl

@@ -155,6 +155,9 @@ end
 # apply operators
 include("apply.jl")
 
+# initialization
+include("initialize.jl")
+
 # Verbosity levels
 const WARN_LEVEL = 1
 const STARTSTOP_LEVEL = 2

src/apply.jl

@@ -1,3 +1,10 @@
+@doc """
+    apply(operator, x)
+
+Apply the operator `operator` to the vector `x`, returning the result.
+By default this falls back to `operator * x` for `AbstractMatrix`, and `f(x)` in other cases.
+""" apply
+
 apply(A::AbstractMatrix, x::AbstractVector) = A * x
 apply(f, x) = f(x)
 
@@ -11,6 +18,22 @@ function apply(operator, x, α₀, α₁)
 end
 
 # GKL, SVD, LSMR
+@doc """
+    apply_normal(operator, x)
+
+Apply the operator to the vector `x`, returning the result.
+By default this falls back to `operator * x` for `AbstractMatrix`, `operator[1](x)` when
+the input is a tuple of functions, and `operator(x, Val(false))` for other cases.
+""" apply_normal
+
+@doc """
+    apply_adjoint(operator, x)
+
+Apply the adjoint of the operator to the vector `x`, returning the result.
+By default this falls back to `operator' * x` for `AbstractMatrix`, `operator[2](x)` when
+the input is a tuple of functions, and `operator(x, Val(true))` for other cases.
+""" apply_adjoint
+
 apply_normal(A::AbstractMatrix, x::AbstractVector) = A * x
 apply_adjoint(A::AbstractMatrix, x::AbstractVector) = A' * x
 apply_normal((f, fadjoint)::Tuple{Any,Any}, x) = f(x)

src/deprecated.jl

@@ -4,3 +4,13 @@ Base.@deprecate(basis(F::GKLFactorization, which::Symbol), basis(F, Val(which)))
 
 import LinearAlgebra: mul!
 Base.@deprecate(mul!(y, b::OrthonormalBasis, x::AbstractVector), unproject!!(y, b, x))
+
+Base.@deprecate(eigsolve(A::AbstractMatrix, howmany::Int, which::Selector, T::Type;
+                         kwargs...),
+                eigsolve(A, Random.rand!(similar(A, T, size(A, 1))), howmany, which;
+                         kwargs...),
+                false)
+
+Base.@deprecate(eigsolve(f, n::Int, howmany::Int, which::Selector, T::Type; kwargs...),
+                eigsolve(f, Random.rand(T, n), howmany, which; kwargs...),
+                false)

src/eigsolve/eigsolve.jl

@@ -12,13 +12,15 @@ the function `f`. Return eigenvalues, eigenvectors and a `ConvergenceInfo` struc
 
 The linear map can be an `AbstractMatrix` (dense or sparse) or a general function or
 callable object. If an `AbstractMatrix` is used, a starting vector `x₀` does not need to be
-provided, it is then chosen as `rand(T, size(A, 1))`. If the linear map is encoded more
-generally as a a callable function or method, the best approach is to provide an explicit
-starting guess `x₀`. Note that `x₀` does not need to be of type `AbstractVector`; any type
-that behaves as a vector and supports the required methods (see KrylovKit docs) is accepted.
-If instead of `x₀` an integer `n` is specified, it is assumed that `x₀` is a regular vector
-and it is initialized to `rand(T, n)`, where the default value of `T` is `Float64`, unless
-specified differently.
+provided, it is then chosen as `rand!(similar(A, T, size(A, 1)))`. If the linear map is
+encoded more generally as a a callable function or method, the best approach is to provide
+an explicit starting guess `x₀`. Alternatively, you can overload [`eigsolve_init`](@ref) to
+construct a suitable starting vector for the function `f`.
+Note that `x₀` does not need to be of type `AbstractVector`; any type that behaves as a
+vector and supports the required methods (see KrylovKit docs) is accepted. If instead of
+`x₀` an integer `n` is specified, it is assumed that `x₀` is a regular vector and it is
+initialized to `rand(T, n)`, where the default value of `T` is `Float64`, unless specified
+differently.
 
 The next arguments are optional, but should typically be specified. `howmany` specifies how
 many eigenvalues should be computed; `which` specifies which eigenvalues should be
@@ -181,20 +183,12 @@ EigSorter(f::F; rev=false) where {F} = EigSorter{F}(f, rev)
 
 const Selector = Union{Symbol,EigSorter}
 
-function eigsolve(A::AbstractMatrix,
-                  howmany::Int=1,
-                  which::Selector=:LM,
-                  T::Type=eltype(A);
-                  kwargs...)
-    x₀ = Random.rand!(similar(A, T, size(A, 1)))
-    return eigsolve(A, x₀, howmany, which; kwargs...)
+function eigsolve(f, howmany::Int, which::Selector=:LM; kwargs...)
+    x₀ = initialize_vector(eigsolve, f)
+    return eigsolve(f, x₀, howmany, which; kwargs...)
 end
-
-function eigsolve(f, n::Int, howmany::Int=1, which::Selector=:LM, T::Type=Float64;
+function eigsolve(f, x₀=initialize_vector(eigsolve, f), howmany::Int=1, which::Selector=:LM;
                   kwargs...)
-    return eigsolve(f, rand(T, n), howmany, which; kwargs...)
-end
-function eigsolve(f, x₀, howmany::Int=1, which::Selector=:LM; kwargs...)
     Tx = typeof(x₀)
     Tfx = Core.Compiler.return_type(apply, Tuple{typeof(f),Tx})
     T = Core.Compiler.return_type(dot, Tuple{Tx,Tfx})

src/initialize.jl

@@ -0,0 +1,18 @@
+"""
+    initialize_vector(f, A)
+
+Construct a starting vector for the Krylov subspace of the operator `A` constructed for `f`.
+For `A::AbstractMatrix`, the default is a random vector of the same size as the number of rows of `A`.
+For a function `A` or custom object, you should either provide a starting vector `x₀` or implement this function.
+"""
+function initialize_vector(f, A::AbstractMatrix)
+    return Random.rand!(similar(A, scalartype(A), size(A, 1)))
+end
+
+function initialize_vector(f, A)
+    error("""
+          Cannot construct a starting vector for the Krylov subspace from the given operator.
+          Either provide a starting vector `x₀`, or implement [`initialize_vector`](@ref) for values of type `$(typeof(f)), $(typeof(A))`.
+          """)
+    return nothing
+end