Block Lanczos for Multiple eigenvalues

Project.toml

@@ -8,6 +8,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 PackageExtensionCompat = "65ce6f38-6b18-4e1d-a461-8949797d7930"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 VectorInterface = "409d34a3-91d5-4945-b6ec-7529ddf182d8"
 
 [weakdeps]
@@ -26,6 +27,7 @@ Logging = "1"
 PackageExtensionCompat = "1"
 Printf = "1"
 Random = "1"
+SparseArrays = "1.11.0"
 Test = "1"
 TestExtras = "0.2,0.3"
 VectorInterface = "0.5"

docs/src/man/Block_Lanczos.md

@@ -0,0 +1,43 @@
+# General Block Lanczos
+$A$ is a Hermitian matrix, 
+$$A = Q T Q'$$
+
+$T$ is a tridiagonal matrix, 
+$$T = \begin{pmatrix}
+    \alpha_1 & \beta_1 \\
+    \beta_1 & \alpha_2 & \beta_2 \\
+    & \beta_2 & \ddots & \ddots \\
+    && \ddots & \alpha_{n-1} & \beta_{n-1} \\
+    &&& \beta_{n-1} & \alpha_n
+\end{pmatrix}$$
+
+But $\beta \neq 0$ and thus eigenvalues of $T$ are different from each other. 
+
+To get multiple eigenvalues, we can use block Lanczos. 
+
+$$A = Q T Q', \quad Q = [X_1|..|X_n],\quad \text{size}(X_i) = (p,p)$$
+
+$$T = \begin{pmatrix}
+    M_1 & B_1'\\
+    B_1 & M_2 & B_2'\\
+    & B_2 & \ddots & \ddots \\
+    && \ddots & M_{n-1} & B_{n-1}'\\
+    &&& B_{n-1} & M_n
+\end{pmatrix}$$
+
+It's advantage is to get multiple eigenvalues. But it's disadvantage is that it can cause ghost eigenvalues because of the loss of orthogonality of the Lanczos vectors. 
+
+So my solution is to make sure each $X_i$ we get is orthogonal to $X_{i_1},..,X_1$:
+
+$$X_i = X_i - Q_{i-1}*(Q_{i-1}'X_i)\\
+Q_{i-1} = [X_1|..|X_{i-1}]
+$$
+
+I use Modified Schidi's method to force the orthogonality of the basis And use some skills to improve the method and speed up. 
+
+# Difference between Block Lanczos and Lanczos in code
+
+1. I don't use inner because it maps a couple of matrix to a scalar. 
+2. I have add test to make sure Block Lanczos can work for map input
+3. I add SaprseArray to do test for sparse matrix input
+

src/algorithms.jl

@@ -108,6 +108,7 @@ Use `Arnoldi` for non-symmetric or non-Hermitian linear operators.
 See also: `factorize`, `eigsolve`, `exponentiate`, `Arnoldi`, `Orthogonalizer`
 """
 struct Lanczos{O<:Orthogonalizer,S<:Real} <: KrylovAlgorithm
+    block_size::Int
     orth::O
     krylovdim::Int
     maxiter::Int
@@ -116,13 +117,14 @@ struct Lanczos{O<:Orthogonalizer,S<:Real} <: KrylovAlgorithm
     verbosity::Int
 end
 function Lanczos(;
+                 block_size::Int=1,
                  krylovdim::Int=KrylovDefaults.krylovdim[],
                  maxiter::Int=KrylovDefaults.maxiter[],
                  tol::Real=KrylovDefaults.tol[],
                  orth::Orthogonalizer=KrylovDefaults.orth,
                  eager::Bool=false,
                  verbosity::Int=KrylovDefaults.verbosity[])
-    return Lanczos(orth, krylovdim, maxiter, tol, eager, verbosity)
+    return Lanczos(block_size, orth, krylovdim, maxiter, tol, eager, verbosity)
 end
 
 """

src/apply.jl

@@ -1,4 +1,4 @@
-apply(A::AbstractMatrix, x::AbstractVector) = A * x
+apply(A::AbstractMatrix, x::AbstractVecOrMat) = A * x
 apply(f, x) = f(x)
 
 function apply(operator, x, α₀, α₁)

src/eigsolve/lanczos.jl

@@ -4,6 +4,9 @@ function eigsolve(A, x₀, howmany::Int, which::Selector, alg::Lanczos;
                                     maxiter=alg.maxiter,
                                     eager=alg.eager,
                                     orth=alg.orth))
+    if alg.block_size > 1
+        return block_lanczos_reortho(A, x₀, howmany, which, alg)
+    end
     krylovdim = alg.krylovdim
     maxiter = alg.maxiter
     if howmany > krylovdim
@@ -150,3 +153,103 @@ function eigsolve(A, x₀, howmany::Int, which::Selector, alg::Lanczos;
            vectors,
            ConvergenceInfo(converged, residuals, normresiduals, numiter, numops)
 end
+
+
+function block_lanczos_reortho(A, x₀::AbstractMatrix{S}, howmany::Int, which::Selector, alg::Lanczos) where S
+    block_size = alg.block_size
+    maxiter = alg.maxiter
+    tol = alg.tol
+    verbosity = alg.verbosity
+    n = size(x₀, 1)
+    max_blocks = n ÷ block_size
+
+    iter = BlockLanczosIterator(A, x₀, block_size, maxiter, alg.orth)
+    fact = initialize(iter; verbosity = verbosity)
+    numops = 2
+
+    converge_check = max(1, 100 ÷ block_size)
+
+    local values, vectors, residuals, normresiduals, num_converged
+    converged = false
+
+	function _res(fact, A, howmany, tol, block_size)
+		T = triblockdiag(fact)
+
+		D, U = eigen(Hermitian((T+T')/2))
+
+		by, rev = eigsort(which)
+		p = sortperm(D; by = by, rev = rev)
+		D = D[p]
+		U = U[:, p]
+
+		howmany_actual = min(howmany, length(D))
+		values = D[1:howmany_actual]
+
+		basis_so_far = view(fact.V, :, 1:fact.k*block_size)
+		vectors = Vector{typeof(fact.V[:, 1])}(undef, howmany_actual)
+
+		for i in 1:howmany_actual
+			vectors[i] = similar(basis_so_far, size(basis_so_far, 1))
+			mul!(vectors[i], basis_so_far, view(U, :, i))
+		end
+
+		residuals = Vector{typeof(vectors[1])}(undef, howmany_actual)
+		normresiduals = Vector{Float64}(undef, howmany_actual)
+
+		for i in 1:howmany_actual
+			residuals[i] = apply(A, vectors[i])
+			axpy!(-values[i], vectors[i], residuals[i])  # residuals[i] -= values[i] * vectors[i]
+			normresiduals[i] = norm(residuals[i])
+		end
+
+		num_converged = count(nr -> nr <= tol, normresiduals)
+		return values, vectors, residuals, normresiduals, num_converged
+	end
+
+    for numiter in 2:min(maxiter, max_blocks - 2)
+        expand!(iter, fact; verbosity = verbosity)
+        numops += 1
+
+        if (numiter % converge_check == 0) || (fact.normR < tol)
+            values, vectors, residuals, normresiduals, num_converged = 
+                _res(fact, A, howmany, tol, block_size)
+
+            if verbosity >= EACHITERATION_LEVEL
+                @info "Block Lanczos eigsolve in iteration $numiter: $num_converged values converged, normres = $(normres2string(normresiduals[1:min(howmany, length(normresiduals))]))"
+            end
+
+			# This convergence condition refers to https://www.netlib.org/utk/people/JackDongarra/etemplates/node251.html
+            if num_converged >= howmany || fact.normR < tol
+                converged = true
+                break
+            end
+        end
+    end
+
+    if !converged
+        values, vectors, residuals, normresiduals, num_converged = 
+            _res(fact, A, howmany, tol, block_size)
+    end
+
+    if (fact.k * block_size > alg.krylovdim)
+        @warn "The real Krylov dimension is $(fact.k * block_size), which is larger than the maximum allowed dimension $(alg.krylovdim)."
+        # In this version we don't shrink the factorization because it might cause issues, different from the ordinary Lanczos.
+        # Why it happens remains to be investigated.
+    end
+
+    if (num_converged < howmany) && verbosity >= WARN_LEVEL
+        @warn """Block Lanczos eigsolve stopped without full convergence after $(fact.k) iterations:
+        * $num_converged eigenvalues converged
+        * norm of residuals = $(normres2string(normresiduals))
+        * number of operations = $numops"""
+    elseif verbosity >= STARTSTOP_LEVEL
+        @info """Block Lanczos eigsolve finished after $(fact.k) iterations:
+        * $num_converged eigenvalues converged
+        * norm of residuals = $(normres2string(normresiduals))
+        * number of operations = $numops"""
+    end
+
+    return values,
+           vectors,
+           ConvergenceInfo(num_converged, residuals, normresiduals, fact.k, numops)
+end