Block Lanczos for Multiple eigenvalues

[yuiyuiui]

I add Bloack Lanczos into eigsove() and add a test for it's accuracy in finding all folds of ground eigenvalues. It shows great stability.

[Jutho]

Thank you for opening this PR. We are definitely interested in having block versions of the different methods available.

However, a general comment about your current implementation is that it is very much tailored to A being some AbstractMatrix subtype and x being a standard Vector. This is very much not the approach that KrylovKit.jl has taken, which tries to be as generic as possible. In particular, vectors can be of arbitrary type, and the only contract is that you can use the methods from VectorInterface.jl, preferably the "potentially mutating" methods with a double bang (!!). VectorInterface.jl provides an implementation of those for standard types from Base, including tuples etc.

The linear map can be an arbitrary function that acts on vectors. In the particular context of block methods, this means for example that they can only act on one vector at the time, i.e. on the different columns of your X matrices.

As a result, the current implementation is really not at all suited to be included in KrylovKit.jl and would need a major rewrite. I can certainly give advise on this, but currently don't have the time to do this myself.

[GiggleLiu]

Thank you for opening this PR. We are definitely interested in having block versions of the different methods available.

However, a general comment about your current implementation is that it is very much tailored to A being some AbstractMatrix subtype and x being a standard Vector. This is very much not the approach that KrylovKit.jl has taken, which tries to be as generic as possible. In particular, vectors can be of arbitrary type, and the only contract is that you can use the methods from VectorInterface.jl, preferably the "potentially mutating" methods with a double bang (!!). VectorInterface.jl provides an implementation of those for standard types from Base, including tuples etc.

The linear map can be an arbitrary function that acts on vectors. In the particular context of block methods, this means for example that they can only act on one vector at the time, i.e. on the different columns of your X matrices.

As a result, the current implementation is really not at all suited to be included in KrylovKit.jl and would need a major rewrite. I can certainly give advise on this, but currently don't have the time to do this myself.

Hi @Jutho , thank you for the suggestion. I can review this PR for the first round. FYI: Implementing block lanczos for KrylovKit is a homework that I assigned to students as a challenge problem in my course. Although this student (of mine) is talented, he is relative new to Julia. I will let you know when this PR is ready for you to review.

[yuiyuiui]

That would be great. I haven't thought it through, so I hope that is not too difficult, as the current eigsolve interface has grown organically and allows –in hindsight– for too many different call signatures. However, specialising on the type of x₀ was never part of the plan as allowing arbitrary vector types was always the goal from the beginning.

If you find a way to make it work, it would be good to already foresee a code path where x₀ isa Block but the map is not known to be hermitian, and then just prints an error that a BlockArnoldi method has not yet been implemented but will be in the future (and I would certainly welcome you doing so if you have the time and interest in doing so).

I am honored to contribute to the implementation of algorithms in KrylovKit.jl. I plan to implement the BlockArnoldi method in July.

In the above commits, I have added support for

eigsolve(f, x₀::Block, [howmany = 1, which = :LM]; kwargs...)

along with the corresponding docstring and tests. For the user’s convenience, I have also added Block(vec) as a constructor for Block{T,S}.

I noticed that Core.Compiler.return_type is used to infer the data type of the inner product. However, this approach may fail for InnerProductVec with a non-trivial dot, i.e., when dot != (x, y) -> x' * y:

julia> using KrylovKit, LinearAlgebra

julia> A = rand(2,2)
2×2 Matrix{Float64}:
 0.68118   0.526311
 0.132412  0.0826394

julia> A = A'*A + I
2×2 Matrix{Float64}:
 1.48154   0.369455
 0.369455  1.28383

julia> ip1 = InnerProductVec(rand(3), (x, y) -> x' * y)
InnerProductVec{var"#1#2", Vector{Float64}}([0.37838962534421194, 0.39745381290616866, 0.22788971208789832], var"#1#2"())

julia> Core.Compiler.return_type(KrylovKit.inner, Tuple{typeof(ip1), typeof(ip1)})
Float64

julia> ip2 = InnerProductVec(rand(3), (x, y) -> x' * A * y)
InnerProductVec{var"#3#4", Vector{Float64}}([0.30882023529528346, 0.5849695825013441, 0.10229414046031216], var"#3#4"())

julia> Core.Compiler.return_type(KrylovKit.inner, Tuple{typeof(ip2), typeof(ip2)})
Any

Therefore, I opted to directly use:

S = typeof(inner(vec[1], vec[1]))

Please let me know if you have any suggestions or if there is a better approach I may have overlooked.

[yuiyuiui]

By the way, I noticed that both issymmetric and ishermitian are used in the code. For example:

https://github.com/Jutho/KrylovKit.jl/blob/master/src/eigsolve/eigsolve.jl#L226-L236)

function eigselector(f,
                     T::Type;
                     issymmetric::Bool=false,
                     ishermitian::Bool=issymmetric && (T <: Real),
                     krylovdim::Int=KrylovDefaults.krylovdim[],
                     maxiter::Int=KrylovDefaults.maxiter[],
                     tol::Real=KrylovDefaults.tol[],
                     orth::Orthogonalizer=KrylovDefaults.orth,
                     eager::Bool=false,
                     verbosity::Int=KrylovDefaults.verbosity[],
                     alg_rrule=nothing)

Since issymmetric is implicitly covered by ishermitian, I was wondering if there is a specific reason for checking both? I tried removing the issymmetric checks, and all tests still passed, but perhaps I’m missing some subtlety.

Just curious and would love to understand the rationale if there's one. Thanks!

[Jutho]

I agree that from a mathematics point of view, only ishermitian ought to be there. The issymmetric keyword is mostly there for users that only work with real matrices/maps and are less familiar with the notion of hermiticity. They might want to specify that their linear map is symmetric, by setting issymmetric=true.

[Jutho]

Thanks for the extensive test suite. I have now reviewed "test/Block.jl" (I think I would prefer a lower case filename, for everything except for the files that define the modules like "src/KrylovKit.jl" and the extensions.

I will still review "test/factorize.jl" and "test/eigsolve.jl" and then I think this is completely finished.

[yuiyuiui]

Thanks for the extensive test suite. I have now reviewed "test/Block.jl" (I think I would prefer a lower case filename, for everything except for the files that define the modules like "src/KrylovKit.jl" and the extensions.

I will still review "test/factorize.jl" and "test/eigsolve.jl" and then I think this is completely finished.

Thank you very much for reviewing my tests so thoroughly. I had very little prior experience with writing tests, and most of mine were written by learning from and referencing your own. Reading your test cases has been extremely instructive for me.

I also agree with your point regarding the use of lowercase for "block.jl". I've now completed the changes according to your latest review suggestions.

[Jutho]

This explicit use of verbosity = 0 and 1 should not happen anymore, but I realise this is literally what is also being done in the normal Lanczos. Maybe you could use the following:

Suggested change

	verbosity = 1
	while fact.k < n
	if verbosity == 1
	@test_logs (:warn,) expand!(iter, fact; verbosity=verbosity)
	verbosity = 0
	else
	@test_logs expand!(iter, fact; verbosity=verbosity)
	verbosity = 1
	end
	end
	verbosity = WARN_LEVEL
	while fact.k < n
	if verbosity == WARN_LEVEL
	@test_logs (:warn,) expand!(iter, fact; verbosity=verbosity)
	verbosity = SILENT_LEVEL
	else
	@test_logs expand!(iter, fact; verbosity=verbosity)
	verbosity = WARN_LEVEL
	end
	end

in combination with defining a new level in "src/KrylovKit.jl" at line 159:

const SILENT_LEVEL = 0

and also updating the Lanczos part of the tests, and replacing line 6 of this test file with

    using KrylovKit: SILENT_LEVEL, WARN_LEVEL, EACHITERATION_LEVEL

[yuiyuiui]

You're absolutely right — this is indeed a more proper and consistent way to write the code. I've gone through the entire KrylovKit.jl codebase, searching for all occurrences of verbosity, and carefully revised each instance to eliminate any direct assignments or comparisons with raw integers.

Now, all uses of verbosity follow a consistent pattern, such as:

verbosity = WARN_LEVEL
verbosity >= WARN_LEVEL

Cases like verbosity = -1 have been replaced with verbosity = SILENT_LEVEL - 1,
and verbosity = 4 with verbosity = EACHITERATION_LEVEL + 1,
while 4+ has been adjusted accordingly to EACHITERATION_LEVEL +.

Let me know if you see any areas that could be further improved!

[yuiyuiui]

In addition, since STARTSTOP_LEVEL is used in the tests, I’ve removed using KrylovKit: from factorize.jl and instead added the following to runtests.jl:

using Test, TestExtras, Logging
using LinearAlgebra, SparseArrays
using KrylovKit
using KrylovKit: SILENT_LEVEL, WARN_LEVEL, STARTSTOP_LEVEL, EACHITERATION_LEVEL
using VectorInterface

Please let me know if there’s a better way to handle this — I’m happy to adjust things as needed.

[Jutho]

Maybe it is useful to have a small utility function that generates some matrix with at least one repeated eigenvalue with a given multiplicity, that can be chosen equal to block_size, something like:

function mat_with_eigrepition(T, N, multiplicity)
    U = qr(randn(T, (N, N))).Q # Haar random matrix
    D = sort(randn(real(T), N))
    i = 0
    while multiplicity >= 2 && (i + multiplicity) <= N
        D[i .+ (1:multiplicity)] .= D[i+1]
        i += multiplicity
        multiplicity -= 1
    end
    A = U * Diagonal(D) * U'
    return (A + A')/2
end

which generates a random Hermitian matrix with eigenvalues with multiplicities multiplicity, multiplicity-1, multiplicity-2, ... in its :SR eigenvalues.

[yuiyuiui]

You're absolutely right — it's best to test BlockLanczos using operators with repeated eigenvalues. I've now replaced the operators in tests of BlockLanczos in both test/eigsolve.jl and test/factorize.jl (except for the one in "Compare Lanczos and BlockLanczos") with ones constructed using mat_with_eigrepition. I didn't replace the operators in "Compare Lanczos and BlockLanczos" because the standard Lanczos method is not capable of computing repeated eigenvalues.

The changes for this round of review have been completed.