Lanczos algorithm + EDKit

[HamidArianZad]

Dear Jie,

Thank you for your excellent package.

I have successfully solved the problem with a spin-1 diamond chain up to 12 spins. Now, I aim to increase the number of spins to 18 while conserving SzT by applying TranslationalBasis. I have 64 GB of RAM, and the code gets terminated when the total Sz exceeds sector 6. To address this issue, I need to utilize the Lanczos algorithm. Here is the code I have managed:

using EDKit, LinearAlgebra
using DelimitedFiles
using KrylovKit
# Hamiltonian parameters:
J1 = 1.0
J2 = J1
J3 = 2.8
J4 = J1
J5 = J1
L = 6
N = 3*L
k=0
sz_min = 0
sz_max = N
B_min = 0
B_max = 0
uB = 0.01
# -------------------------------------------------------
all_energies = []
Size_EigVals = []
for i in sz_min:sz_max
        # onsite Hamiltonian matrix
        O1 = spin((J1, "xx1"), (J1, "yy1"), (J1, "zz1"), D=3)
        O2 = spin((J2, "x1x"), (J2, "y1y"), (J2, "z1z"), D=3)
        O3 = spin((J3, "1xx"), (J3, "1yy"), (J3, "1zz"), D=3)
        # hz = spin((ha, "z11"), (hb, "1z1"), (hc, "11z"), D=3)
        h1 = O1 + O2 + O3
        C1 = spin((J4, "1x1x11"), (J4, "1y1y11"), (J4, "1z1z11"), D=3)
        C2 = spin((J5, "11xx11"), (J5, "11yy11"), (J5, "11zz11"), D=3)
        h2 = C1 + C2
        Nj = diag( spin((1,"z11"), (1,"1z1"), (1,"11z"), D=3) + 3*I )
        total_z(s) = sum(Nj[sj+1] for sj in s)
        basis = TranslationalBasis(L=L, f = s -> total_z(s) == i, k=k, base=3^3)
        H = trans_inv_operator(h1, 1:1, basis) + trans_inv_operator(h2, 1:2, basis)
        if i <= 6
           energy = eigvals(Hermitian(H))
        else
           energy = eigsolve(H, ones(size(H, 1)), howmany=20, which=:SM, T=Float64)
        end
        println("size Eigenvalues = ", size(energy, 1))
        push!(Size_EigVals, size(energy, 1))
        for sublist in energy
            for element in sublist
                push!(all_energies, element)
            end
        end
    println("*************************************************************")
    println("-------------- Finished job: SzT = ", i, "--------------------")
    println("*************************************************************")
end

open("Energies_DiamondS1_cSzT_k0_EDKit_EnDig_J$(J3)_N$(N).txt", "w") do io
    for energies in all_energies
        for En in energies
            writedlm(io, En, '\n')
        end
    end
end

open("Size_EigVals_J$(J3)_N$(N).txt", "w") do io
           writedlm(io, [Size_EigVals], '\t')
end

But I get following error:

LoadError: Tuple field type cannot be Union{}
Stacktrace:
 [1] eigsolve(f::EDKit.Operator{Float64, TranslationalBasis{Int64, Float64}}, x₀::Vector{Float64}, howmany::Int64, which::Symbol; kwargs::@Kwargs{howmany::Int64, which::Symbol, T::DataType})
   @ KrylovKit ~/.julia/packages/KrylovKit/diNbc/src/eigsolve/eigsolve.jl:183
 [2] top-level scope
   @ ~/DiamondS1/EDKit_DiamondS1_FTLM.jl:38
in expression starting at /home/DiamondS1/EDKit_DiamondS1_FTLM.jl:21

It seems that eigsolve is unable to recognize the Hamiltonian matrix, whether it's the full Hamiltonian or the reduced Hamiltonian after applying TranslationalBasis.

I would appreciate your comments and suggestions on this issue.

[jayren3996]

Hi Hamid,

The object Operator currently is not an AbstractMatrix. It can multiply to a vector, though. So for the Lanczos algorithm, the way to go is to use the method

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

Here, the first argument can be the multiplication function x -> H * x.

[HamidArianZad]

Dear Jie,

I performed

x₀ = ones(size(H, 1))
energy = eigsolve(x -> H * x, x₀, howmany=1, which=:LM)

But still get following error

ERROR: LoadError: MethodError: no method matching eigselector(::var"#3#7"{EDKit.Operator{Float64, TranslationalBasis{Int64, Float64}}}, ::Type{Float64}; howmany::Int64, which::Symbol)

Closest candidates are:
  eigselector(::Any, ::Type; issymmetric, ishermitian, krylovdim, maxiter, tol, orth, eager, verbosity) got unsupported keyword arguments "howmany", "which"
   @ KrylovKit ~/.julia/packages/KrylovKit/diNbc/src/eigsolve/eigsolve.jl:205
  eigselector(::AbstractMatrix, ::Type; issymmetric, ishermitian, krylovdim, maxiter, tol, orth, eager, verbosity) got unsupported keyword arguments "howmany", "which"
   @ KrylovKit ~/.julia/packages/KrylovKit/diNbc/src/eigsolve/eigsolve.jl:237

[jayren3996]

Hi Hamid, The "howmany" and "which" are not keyword argument, you can try removing the argument names in the eigsolve function.

[HamidArianZad]

Thank you for your response,

Term energy = eigsolve(x -> Array(H) * x, x₀) still does not work, and I get similar error. Please know that I need to limit the number of eigenvalues to smaller values compared to the hole spectrum specifically for 18 number of spins. So, the term howmany may play important role. I used the alternative term Arpack.eigs :

if i <= 3
        energy = eigvals(Hermitian(H))
    else
        # x₀ = ones(size(H, 1))
        # energy = eigsolve(x -> Array(H) * x, x₀)
        # energy = eigsolve(x -> H * x, x₀)
        energy = Arpack.eigs(Array(H), nev=100, which=:LM, tol=1e-8)

end

Arpack.eigs seems work, but it always returns 6 eigenvalues. I tried other option: which=:SR, but Arpack.eigs does not give more than 6 eigenvalues.

[jayren3996]

Sorry for the confusion. What I mean is

energy = eigsolve(x -> H * x, x₀, 1, :LM)

That is the usage of eigsolve in KrylovKit when the first argument is a function.

Besides, when H is getting large, you can also consider using sparse(H) to get a sparse matrix of the Hamiltonian (and for that purpose you need to import SparseArrays). Using Array(H) is not very efficient in this case.

[HamidArianZad]

I tried what you suggested:

x₀ = ones(size(H, 1))
energy = eigsolve(x -> H * x, x₀, 1, :SR)

still I get error:

ERROR: LoadError: MethodError: no method matching iterate(::KrylovKit.ConvergenceInfo{Vector{Float64}, Vector{Vector{ComplexF64}}})

Closest candidates are:
  iterate(::KrylovKit.SplitRange)
   @ KrylovKit ~/.julia/packages/KrylovKit/diNbc/src/KrylovKit.jl:83
  iterate(::KrylovKit.SplitRange, ::Any)
   @ KrylovKit ~/.julia/packages/KrylovKit/diNbc/src/KrylovKit.jl:83
  iterate(::Core.MethodMatch, ::Int64)
   @ Base deprecated.jl:265

I also considered x = ones(size(H, 1)), but it did not help to solve the problem.

btw, when I consider sparse(H) I get following error from line if i <= 3 energy = eigvals(sparse(H)) :

ERROR: LoadError: MethodError: no method matching eigvals!(::SparseMatrixCSC{Float64, Int64})

Closest candidates are:
  eigvals!(::AbstractMatrix{T}, ::Cholesky{T}; sortby) where T<:Number
   @ LinearAlgebra /snap/julia/92/share/julia/stdlib/v1.10/LinearAlgebra/src/symmetriceigen.jl:217
  eigvals!(::SymTridiagonal{<:Union{Float32, Float64}, <:StridedVector{T} where T})
   @ LinearAlgebra /snap/julia/92/share/julia/stdlib/v1.10/LinearAlgebra/src/tridiag.jl:265
  eigvals!(::SymTridiagonal{<:Union{Float32, Float64}, <:StridedVector{T} where T}, ::UnitRange)
   @ LinearAlgebra /snap/julia/92/share/julia/stdlib/v1.10/LinearAlgebra/src/tridiag.jl:268
  ...

[jayren3996]

Well it is a bit strange. On my computer, for i=7, eigsolve(x -> H * x, x₀, 1, :SR) works fine.

Anyway for Lanczos problem, you can alway convert the Hamiltonian to sparse matrix and then use the Arpack to solve the eigenvalues.

eigs(sparse(H), nev = 20, which=:SR)

As for your second error message. It is just because eigvals function does not work for sparse matrices. So for the complete diagonalization case, you don need to convert H to a sparse matrix.

[HamidArianZad]

Thank you very much for your assistance. The command energy, _ = eigs(sparse(H), nev = 20, which=:SR) works well for system with N <= 12. However, when I attempted to diagonalize the Hamiltonian for the 18-spin system, I encountered a memory error. This happened even when I reduced nev to 1, indicating that the memory issue persists regardless of the number of desired eigenvalues.

The error occurs specifically at i = 7 for the 18-spin system, leading to the termination of the diagonalization process. This is unexpected as I am operating on a PC with 64 gig of RAM. I had the same issue when I tried to use ED method to diagonalize the spin-1 system with N = 18 using ITensor.jl. I am not sure that whether HPhi can do diagonalization for 18-spin system using Lanczos algorithm on my PC.

I would greatly appreciate it if you could provide any insights or suggestions regarding this memory issue. Perhaps there are alternative approaches or optimizations that could be applied by EDKit to address this challenge.