Krylovkit tolerance changes results drastically

I’m writing a difficult to read code, just because I know the analytical result which is -0.274519. This value is meant to be compared with what is printed in the 4th column of this (a laptop is more than enough).

When i change Lanczos tolerance (line 45) i expect something to change. Default for Float64 is 1e-12. I even used 1e-21 and the result still changes on the sixth digit and I’m using Float64.

It seems that the best result (It’s an iterative variational algorithm, so it means it’s converging up to the analytical value) is reached with tol=1e-15. If i put less or more it’s more distant to the value I want it to reach.

Do not mistake what I wrote as tol_zero, that is something that let me supress numerical noise in the matrix and it saturates at 1e-14, which is what i used there.

You can change lanczos tolerance directly from line 45 and see how the result changes.
Lanczos is used on a sparse matrix and i use both eigenalues and eigenvectors.

using Printf
using SparseArrays
using Random
using LinearAlgebra
using BenchmarkTools
using KrylovKit

function IDMRG(tol_zeros::Real, Delta::Real, max_states::Integer, NIter::Integer, Energy::Vector, Energy_2::Vector, Energy_bond::Vector, lattice_size::Vector, Truncation_Error::Vector)


    Sz = sparse(Float64[[0.5, 0] [0, -0.5]])
    Sp = sparse(Float64[[0, 0] [1.0, 0]])
    Sm = sparse(Float64[[0, 1.0] [0, 0]])
    Ham = sparse(Float64[0 0; 0 0])

    Block_Sp = Sp
    Block_Sm = Sm
    Block_Sz = Sz
    Block_Ham = Ham

    m::Int64 = 2 #hilbert size

    temp_matr = Array{Float64}(undef, max_states*2, max_states)

    for i = 1:NIter
        lattice_size[i] = 2*i +2 #spatial row_size      

        #Adding 1 site
        H_new = kron(Block_Ham, sparse(I,2,2)) .- Delta*kron(Block_Sz, Sz) .+ 0.5*(kron(Block_Sp, Sm) .+ kron(Block_Sm, Sp))
        Block_Sz = kron(sparse(I,m,m), Sz) #definition of new Block_Sz
        Block_Sp = kron(sparse(I,m,m), Sp) #definition of new Block_Sp
        Block_Sm = kron(sparse(I,m,m), Sm) #definition of new Block_Sm

        H_super = kron(H_new, sparse(I,m*2,m*2)) .+ kron(sparse(I,m*2,m*2), H_new) .- Delta*kron(Block_Sz, Block_Sz) .+ 0.5*(kron(Block_Sp, Block_Sm) .+ kron(Block_Sm, Block_Sp))
        droptol!(H_super, tol_zeros)

        #trunc index
        m_next::Int64 = m*m*4 # hilbert space row_size SuperBlock
        rho_dim::Int64 = m*2 #sqrt(m_next)
        trunc_index::Int64 = min(rho_dim,max_states) #this is applied to density matrix diag

        #Diagonalizing SuperBlock Lanczos
        Energy_ground, evec_computed, info = eigsolve(H_super, rand(Float64, m_next), 1, :SR, Lanczos(krylovdim=10, tol=1e-15)) # maxiter = 100
        if info.converged < 1
            println("Lanczos did not converge in iteration: ", i)
        Energy[i] = Energy_ground[1] #output of eigensolve is a vector of eigenvals
        psi_ground = evec_computed[1]

        Energy_2[i] = Energy[i]/lattice_size[i]
        if i> 1
            Energy_bond[i] = (Energy[i] - Energy[i-1]) / 2
            Energy_bond[i] = Energy[i]/2

        #Reduced Density Matrix
        psi_matr = reshape(psi_ground, (rho_dim, rho_dim)) #psi matr is a view of psi allocated data
        rho_matr =  psi_matr * psi_matr' #a way of tracing rho_matrix on the right part
        rho_matr_sym = Symmetric(rho_matr) #it seems this is indeed symmetric
        #Diagonalizing Density Matrix, Trunc matrix
        evals, evecs = eigen!(rho_matr_sym, sortby=-) #if I use eigen! i need to not use rho_matrx again

        trunc_evals = @view evals[1:trunc_index]
        trunc_matrx = @view evecs[:,1:trunc_index]
        trunc_matrx .= ifelse.(abs.(trunc_matrx) .< tol_zeros, 0, trunc_matrx)
        Truncation_Error[i] = 1 - sum(trunc_evals)

        # Preallocation for mul!
        temp = @view temp_matr[1:rho_dim,1:trunc_index]

        #Higher dimension kron
        mul!(temp, H_new, trunc_matrx)
        Block_Ham = trunc_matrx' * temp
        Block_Ham .= ifelse.(abs.(Block_Ham) .< tol_zeros, 0, Block_Ham)

        mul!(temp, Block_Sz, trunc_matrx)
        Block_Sz = trunc_matrx' * temp
        Block_Sz .= ifelse.(abs.(Block_Sz) .< tol_zeros, 0, Block_Sz)

        mul!(temp, Block_Sp, trunc_matrx)
        Block_Sp = trunc_matrx' * temp
        Block_Sp .= ifelse.(abs.(Block_Sp) .< tol_zeros, 0, Block_Sp)

        mul!(temp, Block_Sm, trunc_matrx)
        Block_Sm = trunc_matrx' * temp
        Block_Sm .= ifelse.(abs.(Block_Sm) .< tol_zeros, 0, Block_Sm)

        m = trunc_index # hilbert space row_size before next kron products

#Definition of starting parameters
tol_zeros = 1e-14
Delta::Float64 = 0.5
max_states::Int16 = 5
NIter::Int32 = 30 #Final lattice size: 2*NIter + 2

#output vectors
lattice_size = Vector{Int64}(undef, NIter)
Energy = Vector{Float64}(undef, NIter)
Energy_2 = Vector{Float64}(undef, NIter)
Energy_bond = Vector{Float64}(undef, NIter)
Truncation_Error = Vector{Float64}(undef, NIter)

IDMRG(tol_zeros, Delta, max_states, NIter, Energy, Energy_2, Energy_bond, lattice_size, Truncation_Error)

for i = 1:NIter
    print(lpad(lattice_size[i], 2, "0"), " ")
    @printf("| %.15e | %.15e | %.15e | %.15e | \n", Energy[i], Energy_bond[i], Energy_2[i], Truncation_Error[i])