KrylovKit, syntax issue

General Usage
,
1

So, I’m trying to use Lanczos diagnalization of a sparse marix of float64. Now I don’t want to specify the hint x0, so i just use rand():

Energy[i], psi_ground, info = eigsolve(H_super, rand(Float64, m*2), 1, :SR, Lanczos)

Where m*2 is the size of x0.

I’m following the expert notation of Eigenvalue problems · KrylovKit.jl eigensolve.

The error is:

ERROR: LoadError: MethodError: no method matching eigsolve(::Matrix{Float64}, ::Vector{Float64}, ::Int64, ::Symbol, ::Type{Lanczos})

Closest candidates are:
 eigsolve(::Any, ::Int64, ::Int64, ::Union{Symbol, EigSorter}, ::Type; kwargs...)
  @ KrylovKit ~/.julia/packages/KrylovKit/r8GLV/src/eigsolve/eigsolve.jl:178
 eigsolve(::Any, ::Any, ::Int64, ::Union{Symbol, EigSorter}, ::Arnoldi)
  @ KrylovKit ~/.julia/packages/KrylovKit/r8GLV/src/eigsolve/arnoldi.jl:130
 eigsolve(::Any, ::Any, ::Int64, ::Union{Symbol, EigSorter}, ::Lanczos)
  @ KrylovKit ~/.julia/packages/KrylovKit/r8GLV/src/eigsolve/lanczos.jl:1
 ...

P.S: is it really recognizing it as a sparse matrix? I used Kronecker.jl to build the sparse matrix H_super, which doesn’t follow sparseArrays notation.

Another question: is it possible to use kwargs on expert notation? I just want to specify Lanczos algorithm along all the other features

2

Your last argument should be Lanczos(). Or Lanczos(; kwargs...) if you want to specify additional options such as krylovdim, tol and maxiter.

The error signature also indicates that your H_super is a dense matrix of standard type Matrix{Float64}, so if you intended for this to be a sparse array, you did something wrong in constructing it.

3

kronecker from Kronecker.jl is never recognizes by SparseArrays.jl library as sparse. Does it mean that only kron works with KrylovKit?

To make an example I wrote something like:
H_super = kronecker(dense_matrix, sparse(I, 4,4))
Where dense_matrix is a 4*4 matrix

4

Another issue.
I tried using eigsolve to compute the smallest eigenvector, but it seems that even specifying howmany=1, it computes 5 different eigenvectors. I compared it with eigvecs, even if the eigenvectors are represented in a different way.

using LinearAlgebra
using KrylovKit

A = Float64[1 0 0 0 0; 0 1 2 0 0; 1 1 0 0 0; 0 1 3 2 0; 2 0 5 2 1]

V = eigvecs(A)

show(V)
println("\n")

val, vec, info = eigsolve(A, rand(Float64, 5), 1, :SR, Lanczos())

show(vec)
println("\n")
5

This is expected behaviour, the idea being that during the algorithm, there might have been more converged eigenvectors than the ones you requested. In that case, rather than discarding this information, KrylovKit just returns all the converged eigenvalue/vector pairs anyways. Note that this is also mentioned in the docstring of eigsolve:

## Return values:

The return value is always of the form `vals, vecs, info = eigsolve(...)` with

    vals: a Vector containing the eigenvalues, of length at least howmany, but could be longer if more eigenvalues were converged at the same cost. Eigenvalues will be real if Lanczos was used and complex if Arnoldi was used (see below).

    vecs: a Vector of corresponding eigenvectors, of the same length as vals. Note that eigenvectors are not returned as a matrix, as the linear map could act on any custom Julia type with vector like behavior, i.e. the elements of the list vecs are objects that are typically similar to the starting guess x₀, up to a possibly different eltype. In particular for a general matrix (i.e. with Arnoldi) the eigenvectors are generally complex and are therefore always returned in a complex number format. When the linear map is a simple AbstractMatrix, vecs will be Vector{Vector{<:Number}}.

If you really only want the smallest eigenvalue, you can thus just get the first entry of val and vec.

1 Like
6

Indeed, this is intended behaviour. The initial motivation for this was, if you request for example the largest magnitude (LM) eigenvalue, there might be several eigenvalues of equal magnitude. A typical case is a real matrix, where the LM eigenvalue is complex, meaning that actually both the value and its conjugate (real ± imaginary) will be eigenvalues, and have the same magnitude. Such eigenvalues will converge simultaneously, and it wouldn’t be clear which one to return. Hence, instead, it was chosen to return all eigenvalues that have converged at the iteration where the requested number of eigenvalues have converged.

On the other hand, if eigsolve did not manage to converge the requested number of eigenvalues before spending all iterations, it will not return less eigenvalues, but will still return the requested number. They might still be useful, but might just not yet have reached the required precision. With the default verbosity settings, no warning is printed, and the user should check info.converged or normres to see if all eigenvalues have converged to sufficient accuracy.

1 Like
7

what are the default maxiter and krylovdim?
Can’t find anything in the documentation

8

You can find them in the documentation here: Available algorithms · KrylovKit.jl

Alternatively, this is the docstring of KrylovKit.KrylovDefaults, so if you type that in your console in helpmode that should pull up the same information:

julia> using KrylovKit

help?> KrylovKit.KrylovDefaults
  module KrylovDefaults
      const orth = KrylovKit.ModifiedGramSchmidtIR()
      const krylovdim = 30
      const maxiter = 100
      const tol = 1e-12
  end

  A module listing the default values for the typical parameters in Krylov
  based algorithms:

    •  orth = ModifiedGramSchmidtIR(): the orthogonalization routine used
       to orthogonalize the Krylov basis in the Lanczos or Arnoldi
       iteration

    •  krylovdim = 30: the maximal dimension of the Krylov subspace that
       will be constructed

    •  maxiter = 100: the maximal number of outer iterations, i.e. the
       maximum number of times the Krylov subspace may be rebuilt

    •  tol = 1e-12: the tolerance to which the problem must be solved,
       based on a suitable error measure, e.g. the norm of some residual.

  │ Warning
  │
  │  The default value of tol is a Float64 value, if you solve problems
  │  in Float32 or ComplexF32 arithmetic, you should always specify a
  │  new tol as the default value will not be attainable.
9

wait, how come 1e-12 is default for float64? sometimes i have values which reach the 16th decimal digit with Float64

10

This is a consequence of how floating point values work in general.
Note that floating point numbers cannot represent all numbers. In particular, you can ask for the next floating point value, which is not something you can do for the regular real numbers. Usually, the smallest number you can add to a floating point number in order to obtain a number that is not equal, is denoted by eps. (also note that eps is relative to the number you’re considering) For example, the following statements hold:

julia> 1.0 + eps(1.0) == 1.0
false

julia> 1.0 + eps(1.0) / 2 == 1.0
true

For Float64 the values of eps(1.0) is around 1e-16. The error measure of KrylovKit is based around this, taking into account that it is comparing two vectors that are normalized to unity, such that at the very best n * eps() can be expected, where n is the size of the vector, and that there are some machine precision factors that make it so that the best precision is not attainable.

If you want to learn more about this behaviour, the wikipedia page is already quite extensive: Floating-point arithmetic - Wikipedia

11

I don’t know if I should paste my program here, as it doesn’t give the expected results. However it gives results and changing tol changes them. I expect a saturation of precision by changing the tolerance. So if you write tol=1e-22, the results should be the same as for tol=1e-16 for a Float64. I should add that even increasing tol to 1e-19 changes them

12

Have you checked info.converged wether KrylovKit.jl actually thinks your vectors are converged?
Maybe you need to increase krylovdim? Sometimes that is necessary for convergence.

13

I ask just lowest energy. I set krylovdim as 10. The converged are always >=1 (printed info.converged)

14

Ok what do you mean by

?
Changes dramatically? That shouldn’t be. Small changes are of course expected due to finite numerical precision. I wouldn’t expect some saturation. If you request too low tolerance I would expect KrylovKit to say that the results did not converge.

Is your groundstate degenerate perhaps? In that case you’d get some random linear combination of all groundstates and that would change from run to run.

15

it should not. I’m still trying to find out the bug that gives me wrong results (and it may be related to this bug too). However i just encountered issue in mkl, as it changes the results. I used a simpler program without krilovkit to address it.

16

Following the issue in MKL package wrong results - #3 by Enlil50 it seems that wrong result were brought upon little (non)zeros of order 10^-16 roaming around the matrix. Whereas for dense matrices I know 2 methods for ensuring symmetricity:

  1. A_sym = 0.5*(A + A’)
  2. A_sym = Symmetric(A) #this is just a view of the upper triangular section

For sparse matrix I can use just the first one. It seems that it doesn’t change much. just to approach the issue in a different way: is there something which suppresses numbers under a certain tolerance or threshold in a given sparse matrix ,matrix or array? I would like to try it first.

P.S: i expect that in sparse matrix the zeros are treated just like sparse matrix zeros: they should not be considered in linear algebra operations and such. I don’t need to reallocate the vectors, but it could even give me a little boost.

17

Try droptol!(A, tol).

18

It solved the issue. Now i need to know what is the real tolerance I need to use.

Regarding this thread issue: even with and without droptol, they change a lot. I ranged from 1e-12 to 1e-19 and the significant digit which changes is even before the last Float32 one. However I’m using Float64 everywhere.

This is the program, be sure to comment MKL and MKLSparse as they don’t really work. (Lanzos is at line 43 if you copy the program and tol is at line 35)

using Printf
using SparseArrays
using Random
using LinearAlgebra
using BenchmarkTools
# using MKL
# using MKLSparse
using KrylovKit

function IDMRG(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(zeros(2,2))

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

    m::Int64 = 2 #hilbert size

    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

        #SuperBlock
        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, 1e-15)

        #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-16)) # maxiter = 100
        if info.converged < 1
            println("Lanczos did not converge in iteration: ", i)
            exit()
        end
        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
        else
            Energy_bond[i] = Energy[i]/2
        end


        #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]
        Truncation_Error[i] = 1 - sum(trunc_evals)


        # Preallocation for mul!
        temp = Array{Float64}(undef, rho_dim, trunc_index)

        #Higher dimension kron
        mul!(temp, H_new, trunc_matrx)
        Block_Ham = trunc_matrx' * temp

        mul!(temp, Block_Sz, trunc_matrx)
        Block_Sz = trunc_matrx' * temp

        mul!(temp, Block_Sp, trunc_matrx)
        Block_Sp = trunc_matrx' * temp

        mul!(temp, Block_Sm, trunc_matrx)
        Block_Sm = trunc_matrx' * temp


        m = trunc_index # hilbert space row_size before next kron products
    end        
end

Random.seed!(1)

#Definition of starting parameters
Delta::Float64 = 0.5
max_states::Int16 = 5
NIter::Int32 = 15 #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(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("| %.16f | %.16f | %.16f | %.16f | \n", Energy[i], Energy_bond[i], Energy_2[i], Truncation_Error[i])
end
19

Is there something like this for dense matrices too?

20

Not specifically. You can do something like:

A = rand(10,10) # some matrix
tol = 1e-2 # some tolerance
A .= ifelse.(A .< tol, 0, A)
1 Like