Block Lanczos for Multiple eigenvalues #125

yuiyuiui · 2025-03-24T16:02:22Z

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 · 2025-03-28T09:54:49Z

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 · 2025-03-29T09:24:45Z

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.

Block_Lanczos.md

Project.toml

src/algorithms.jl

test/runtests.jl

test/eigsolve.jl

src/algorithms.jl

Project.toml

docs/src/man/Block_Lanczos.md

…read

…t space)

src/factorizations/blocklanczos.jl

test/Block.jl

test/ad/eigsolve.jl

Jutho · 2025-05-27T10:15:22Z

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 · 2025-05-27T14:37:01Z

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.

test/eigsolve.jl

Jutho · 2025-05-27T23:25:48Z

test/factorize.jl

+            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


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

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!

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 · 2025-05-27T23:37:44Z

test/factorize.jl

+    using KrylovKit: EACHITERATION_LEVEL
+    @testset for T in scalartypes
+        A = rand(T, (N, N))
+        A = (A + A') / 2


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.

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.

yuiyuiui · 2025-06-06T15:45:48Z

@Jutho Thanks for your previous review, Professor Jutho. I've implemented all the suggested changes.

Also, during a recent check, I noticed that I had mistakenly referred to the removed BlockKrylovFactorization type in the docstring — this has now been corrected to KrylovFactorization.

"""
    mutable struct BlockLanczosFactorization{T,S<:Number,SR<:Real} <: KrylovFactorization{T,S,SR}

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

```julia
A * V = V * H + R * B'

Jutho · 2025-06-07T22:56:31Z

test/eigsolve.jl

+    scalartypes = mode === :vector ? (Float32, Float64, ComplexF32, ComplexF64) :
+                  (ComplexF64,)
+    @testset for T in scalartypes
+        block_size = 2


It might be useful to test block_size = 4 for the "large-scale"/iterative example. Then the matrix A constructed will have an eigenvalue with multiplicity 4, one with multiplicity 3, one with multiplicity 2, and all the rest probably multiplicity 1. Since n=10, asking for n eigenvalues, wil include all eigenvalues with multiplicity > 1, at least for the :SR case. Maybe the mat_with_eigrepition can be generalised to also include eigenvalues with multiplicities at the :LR end of the spectrum.

Thank you for the suggestion. I’ve amended the mat_with_eigrepition function so that it can now include eigenvalues with specified multiplicities at the :LR end of the spectrum:

julia> A = mat_with_eigrepition(Float64, 20, 4) 20×20 Matrix{Float64}: 1.64765 0.23796 0.190415 0.0450364 0.0406828 -0.288017 -0.0258486 … -0.055608 -0.0550633 0.248672 -0.0968521 -0.230212 0.0789045 -0.00930553 0.23796 -0.105748 -0.398509 0.516642 -0.101961 0.413682 -0.29606 0.0315604 -0.41862 -0.18788 0.331702 -0.00655033 -0.392029 -0.102039 0.190415 -0.398509 -0.0866956 0.220127 -0.287181 0.283803 0.26477 -0.390588 -0.642699 -0.13838 -0.194419 -0.0671712 0.18917 -0.465996 0.0450364 0.516642 0.220127 0.856838 -0.269163 -0.315103 -0.0386052 -0.11614 0.0742843 -0.319176 0.0529575 0.253753 -0.0365612 -0.0667203 0.0406828 -0.101961 -0.287181 -0.269163 0.594282 0.00163995 0.170808 0.207082 -0.1562 -0.19428 0.415778 -0.368435 -0.230373 -0.171111 -0.288017 0.413682 0.283803 -0.315103 0.00163995 0.0074906 -0.0713997 … 0.12658 0.0077557 0.103292 -0.00558347 -0.34428 0.0530062 -0.468729 -0.0258486 -0.29606 0.26477 -0.0386052 0.170808 -0.0713997 0.159957 -0.151803 0.0611415 9.49214e-5 0.485483 0.634188 0.0777522 -0.0298986 -0.0920549 0.479077 0.0183663 0.326614 -0.0655997 -0.106527 -0.109865 -0.28146 0.146958 0.0442693 -0.0695585 -0.58645 0.754562 0.062727 0.0133804 -0.0585919 0.157328 -0.0633693 -0.346974 -0.159705 -0.325247 -0.00909591 0.302009 -0.0396493 -0.161184 0.0116409 0.48292 -0.52907 -0.344266 -0.0438307 0.0323836 0.166902 -0.181455 -0.193587 -0.0193103 -0.528469 -0.0882282 -0.0320168 0.257896 -0.249788 -0.0226747 -0.0438843 -0.580531 0.408865 0.178618 0.335627 -0.122416 0.00285734 0.315202 … 0.120869 -0.114879 -0.408442 -0.337987 0.259911 -0.349919 -0.526433 0.0583847 0.17269 -0.761463 0.328175 -0.313859 -0.558264 -0.303743 -0.170437 0.158015 0.207404 0.0214014 -0.486985 0.264334 0.130286 0.184077 -0.231295 0.0304647 -0.0288624 -0.193352 0.0903255 0.160203 -0.619488 -0.449529 0.329883 -0.607205 -0.266004 0.00742897 0.144336 -0.055608 0.0315604 -0.390588 -0.11614 0.207082 0.12658 -0.151803 0.674946 0.0106471 -0.00133066 -0.330256 -0.219263 -0.228415 -0.0186149 -0.0550633 -0.41862 -0.642699 0.0742843 -0.1562 0.0077557 0.0611415 0.0106471 0.770097 -0.274651 -0.00239458 0.160631 0.290718 0.240263 0.248672 -0.18788 -0.13838 -0.319176 -0.19428 0.103292 9.49214e-5 … -0.00133066 -0.274651 -0.0845077 -0.204312 -0.329219 -0.315231 -0.304113 -0.0968521 0.331702 -0.194419 0.0529575 0.415778 -0.00558347 0.485483 -0.330256 -0.00239458 -0.204312 0.816954 -0.432552 0.00884376 0.267609 -0.230212 -0.00655033 -0.0671712 0.253753 -0.368435 -0.34428 0.634188 -0.219263 0.160631 -0.329219 -0.432552 0.628372 0.199302 -0.379576 0.0789045 -0.392029 0.18917 -0.0365612 -0.230373 0.0530062 0.0777522 -0.228415 0.290718 -0.315231 0.00884376 0.199302 0.233437 -0.370168 -0.00930553 -0.102039 -0.465996 -0.0667203 -0.171111 -0.468729 -0.0298986 -0.0186149 0.240263 -0.304113 0.267609 -0.379576 -0.370168 -0.171266 julia> eigen(A).values 20-element Vector{Float64}: -1.5470214795656485 -1.5470214795656436 -1.54702147956564 -1.547021479565638 -0.16956006531619608 -0.1695600653161953 -0.16956006531619444 0.028156699450318217 0.02815669945032057 0.2535974736449158 0.345679848482243 0.43944072498824177 0.43944072498824377 1.5804212584185884 1.5804212584185886 1.5804212584185906 2.2669266520772786 2.26692665207728 2.2669266520772804 2.266926652077282

And I’ve also updated the block_size to 4 in the testset "BlockLanczos – eigsolve iteratively".

@testset "BlockLanczos - eigsolve iteratively $mode" for mode in (:vector, :inplace, :outplace) scalartypes = mode === :vector ? (Float32, Float64, ComplexF32, ComplexF64) : (ComplexF64,) @testset for T in scalartypes block_size = 4

test/eigsolve.jl

Jutho · 2025-06-07T23:08:03Z

test/eigsolve.jl

+        evals1, _, info1 = eigsolve(wrapop(A, Val(mode)), x₀_lanczos, n, :SR, alg1)
+        evals2, _, info2 = eigsolve(wrapop(A, Val(mode)), x₀_block, n, :SR, alg2)
+        @test info1.converged >= info2.converged + 1
+    end


Are these two code blocks the same except for block_size? It would probably make sense to just have a for block_size in (1, 4) loop around a single block then?

The code blocks for the two different block_size values are not exactly the same. For block_size = 1, I wanted to demonstrate that BlockLanczos behaves the same as Lanczos, whereas for block_size > 1, BlockLanczos tends to converge more slowly.

That said, I realize this intention wasn't clearly reflected in the code, so I’ve added some additional tests in the block_size = 1 block to make this equivalence between BlockLanczos and Lanczos more explicit.

@testset "Compare Lanczos and BlockLanczos $mode" for mode in (:vector, :inplace, :outplace) scalartypes = mode === :vector ? (Float32, Float64, ComplexF32, ComplexF64) : (ComplexF64,) mode = :inplace T = ComplexF64 @testset for T in scalartypes A = rand(T, (2N, 2N)) .- one(T) / 2 A = (A + A') / 2 block_size = 1 x₀_block = Block{T}([wrapvec(rand(T, 2N), Val(mode)) for _ in 1:block_size]) x₀_lanczos = x₀_block[1] alg1 = Lanczos(; krylovdim=2n, maxiter=10, tol=tolerance(T), verbosity=WARN_LEVEL) alg2 = BlockLanczos(; krylovdim=2n, maxiter=10, tol=tolerance(T), verbosity=WARN_LEVEL) evals1, vectors1, info1 = eigsolve(wrapop(A, Val(mode)), x₀_lanczos, n, :SR, alg1) evals2, vectors2, info2 = eigsolve(wrapop(A, Val(mode)), x₀_block, n, :SR, alg2) @test info1.converged == info2.converged @test info1.numiter == info2.numiter @test info1.numops == info2.numops @test isapprox(info1.normres, info2.normres, atol=tolerance(T)) @test isapprox(unwrapvec.(info1.residual[1:info1.converged])[2], unwrapvec.(info2.residual[1:info2.converged])[2]; atol=tolerance(T)) @test isapprox(evals1[1:info1.converged], evals2[1:info2.converged]; atol= tolerance(T)) end @testset for T in scalartypes A = rand(T, (2N, 2N)) .- one(T) / 2 A = (A + A') / 2 block_size = 4 x₀_block = Block{T}([wrapvec(rand(T, 2N), Val(mode)) for _ in 1:block_size]) x₀_lanczos = x₀_block[1] alg1 = Lanczos(; krylovdim=2n, maxiter=10, tol=tolerance(T), verbosity=WARN_LEVEL) alg2 = BlockLanczos(; krylovdim=2n, maxiter=10, tol=tolerance(T), verbosity=WARN_LEVEL) evals1, _, info1 = eigsolve(wrapop(A, Val(mode)), x₀_lanczos, n, :SR, alg1) evals2, _, info2 = eigsolve(wrapop(A, Val(mode)), x₀_block, n, :SR, alg2) @test info1.converged >= info2.converged + 1 end end

Jutho · 2025-06-07T23:12:19Z

Thanks for the reminder, and for fixing all the verbosity. I did not realize how many instances there were, and in retrospect would probably not have asked it because it also clutters up the PR, but I am very grateful that you took the time to do this and so I suggest to now keep this in this PR.

I left two small final comments for the tests in test/eigsolve.jl. When these are addressed, I think this is all set for being merged.

Jutho · 2025-06-08T22:06:55Z

This is ready to be merged. Thanks for this massive amount of great work and the many useful interactions. This constitutes a great start for more block methods in KrylovKit.jl.

Jutho · 2025-07-31T10:55:21Z

@yuiyuiui , I am currently also revising the BlockLanczos tests a little bit. In the testset on "Compare Lanczos and BlockLanczos", the verbosity level was still set to WARN_LEVEL, which I don't think is necessary as it prints out a lot of warning in the CI process.

I had a question about the comparison where you use BlockLanczos with a block size of 4, so the second part of this test set. In these tests, with those configurations of krylovdim and maxiter, it seems that Lanczos often manages to converge a few eigenvalues (fewer than the requested number, hence a warning), but typically BlockLanczos doesn't manage to get any eigenvalue converged. The only test you have in this case is

@test info1.converged >= info2.converged + 1

with info1 being Lanczos and info2 being BlockLanczos. Why do you expect that Lanczos will always converge at least one more? To me, this is not obvious.

yuiyuiui · 2025-07-31T13:40:32Z

@Jutho I set this test to prove that when block size is more than 1, Blcok Lanczos converge more slowly than Lanczos. The Mathematical theorem refers to Golub, Gene H. “Matrix Computation.”: Corollary 10.1.3 and Theorem 10.3.2.

Take solving the biggest eigenvalue (in norm) as an example. After $k$ iterations, let the largest eigenvalues (in norm) obtained from the original matrix, Lanczos, and Block Lanczos be $\lambda_1$, $\theta_1$, and $\mu_1$, respectively. Then

$$e = \frac{|\lambda_1-\theta_1|}{\lambda_1-\lambda_n} = O( \frac{1}{c_{k-1}(1+2\rho_1)} ),~~ \rho_1 = \frac{\lambda_1-\lambda_2}{\lambda_2-\lambda_n}$$

$$e' = \frac{|\lambda_1-\mu_k|}{\lambda_1-\lambda_n} = O( \frac{1}{c_{k-1}(1+2\rho_1')} ), ~~\rho_1' = \frac{\lambda_1-\lambda_{p+1}}{\lambda_{p+1}-\lambda_n}$$

$c_k$ is the Chebyshev polynomials. When $p>1$, I want to show that, with the same number of Ritz eigenvalues, $e<e'$.

First:

$$\rho_1 = \frac{\lambda_1-\lambda_2}{\lambda_2-\lambda_n}\geq\rho_1' = \frac{\lambda_1-\lambda_{p+1}}{\lambda_{p+1}-\lambda_n}$$

Second: the number of iterations of Lanczos $k$ is greater than $k'$ of Block Lanczos. $k\simeq pk'$

Finaly: when $x>1$, $c_k(x)$ is increasing with $x$ and $c_k'(x)$ increase fast with $k$. So:

$$c_k(1+2\rho_1)>c_{k'}(1+2\rho_1)\geq c_{k'}(1+2\rho_1')$$

$$\Longrightarrow: e < e' $$

I want to show that in test. May I should add a commet in the test: Larger block size leads to slower convergence.

yuiyuiui marked this pull request as draft March 24, 2025 16:20

yuiyuiui changed the title ~~Block Lanczos for Multiple eigenvalues~~ [WIP] Block Lanczos for Multiple eigenvalues Mar 24, 2025

yuiyuiui marked this pull request as ready for review March 24, 2025 16:25

GiggleLiu reviewed Mar 30, 2025

View reviewed changes

yuiyuiui closed this Mar 30, 2025

yuiyuiui force-pushed the master branch from b7e48ca to 68ced46 Compare March 30, 2025 18:33

I have remove all tab

c734c45

yuiyuiui reopened this Mar 30, 2025

GiggleLiu reviewed Mar 31, 2025

View reviewed changes

Project.toml Outdated Show resolved Hide resolved

docs/src/man/Block_Lanczos.md Outdated Show resolved Hide resolved

yuiyuiui and others added 19 commits March 31, 2025 14:23

remove blocklanczos.md, all tabs and polish codes to be more easy to …

adbaf1a

…read

review

7ad9583

inner product version

055c154

mutable block lanczos and tests to be added

29447c7

save

bf0b9fa

add test for abstract_qr

3292f8d

the last test of block lanczos facorizations to do

040577d

format

6156985

revise format

e4c2d64

format revise

6ea0000

save

0a10ca6

format

a2febc6

format

fbaf863

revise project.toml

d5a4618

I am ready

ac130a2

firmat

bb7ea72

add interface for single vector (element in an abstract inner produnc…

a9f0d13

…t space)

revise Project.toml

72d7dbd

add references

706f7c8