Feature request: block lanczos for solving degenerate ground states

[GiggleLiu]

I used the following code for demonstrating the 4-fold degeneracy of toric code model, without success. I leave an issue here so that someone (including me) could inspect it in the future.

using Yao

function toric_code_strings(m::Int, n::Int)
    li = LinearIndices((m, n))
    bottom(i, j) = li[mod1(i, m), mod1(j, n)] + m * n
    right(i, j) = li[mod1(i, m), mod1(j, n)]
    xstrings = Vector{Int}[]
    zstrings = Vector{Int}[]
    for i=1:m, j=1:n
        # face center
        push!(xstrings, [bottom(i, j-1), right(i, j), bottom(i, j), right(i-1, j)])
        # cross
        push!(zstrings, [right(i, j), bottom(i, j), right(i, j+1), bottom(i+1, j)])
    end
    return xstrings, zstrings
end

function toric_code_hamiltonian(m::Int, n::Int)
    xstrings, zstrings = toric_code_strings(m, n)
    sum([kron(2m*n, [i=>X for i in xs]...) for xs in xstrings[1:end-1]]) + sum([kron(2m*n, [i=>Z for i in zs]...) for zs in zstrings[1:end-1]])
end

h = toric_code_hamiltonian(3, 3)

using KrylovKit
eigsolve(-mat(h), 10, :SR)

I got the following output:

julia> eigsolve(-mat(h), 10, :SR)[1]
17-element Vector{Float64}:
 -15.999999999999964
 -13.999999999999982
 -11.999999999999954
  -9.999999999999996
  -8.000000000000004
  -5.999999999999954
  -3.9999999999999893
  -2.0000000000000675
  -2.1316282072803006e-14
   2.0000000000000604
   4.000000000000041
   5.999999999999986
   7.999999999999975
   9.999999999999996
  11.999999999999982
  14.000000000000005
  16.000000000000007

I hope block Lanczos could be supported in the future so that degeneracies could be resolved correctly. For someone who might be interested in implementing this, please check:

• Matrix Computation, by Golub, Sec. 10.3.6

[tibidabo]

Hello why is issue this closed? I would also be interested in a Block Lanczos version.

[Jutho]

It is not closed; the issue FFS-julia#28 was closed.

[hz-xiaxz]

Hello, for your case, I tried ArnoldiMethod.jl and get the output

julia> using ArnoldiMethod
julia> decomp, history  = partialschur(-mat(h),nev=10,which=:SR)
julia> decomp.eigenvalues
10-element Vector{ComplexF64}:
 -16.000000000000053 + 2.6842207305417474e-16im
  -16.00000000000003 - 6.6951290300828266e-18im
 -14.000000000000016 + 1.3094008837548432e-15im
 -14.000000000000012 + 2.78841710585368e-17im
 -12.000000000000014 - 2.796888294211598e-17im
 -12.000000000000005 + 7.040218273262716e-16im
 -10.000000000000005 + 2.315044958024994e-16im
  -9.999999999999998 + 1.5382714567093696e-17im
  -8.000000000000004 - 2.1014552699603364e-16im
  -8.000000000000002 + 1.5056230642491578e-17im

Is that what you want? @GiggleLiu

[GiggleLiu]

Hello, for your case, I tried ArnoldiMethod.jl and get the output

julia> using ArnoldiMethod
julia> decomp, history  = partialschur(-mat(h),nev=10,which=:SR)
julia> decomp.eigenvalues
10-element Vector{ComplexF64}:
 -16.000000000000053 + 2.6842207305417474e-16im
  -16.00000000000003 - 6.6951290300828266e-18im
 -14.000000000000016 + 1.3094008837548432e-15im
 -14.000000000000012 + 2.78841710585368e-17im
 -12.000000000000014 - 2.796888294211598e-17im
 -12.000000000000005 + 7.040218273262716e-16im
 -10.000000000000005 + 2.315044958024994e-16im
  -9.999999999999998 + 1.5382714567093696e-17im
  -8.000000000000004 - 2.1014552699603364e-16im
  -8.000000000000002 + 1.5056230642491578e-17im

Is that what you want? @GiggleLiu

It seems not all minimum eigenvalues are resolved. There should be a 4 fold degeneracy.

[hz-xiaxz]

Hello, for your case, I tried ArnoldiMethod.jl and get the output

julia> using ArnoldiMethod
julia> decomp, history  = partialschur(-mat(h),nev=10,which=:SR)
julia> decomp.eigenvalues
10-element Vector{ComplexF64}:
 -16.000000000000053 + 2.6842207305417474e-16im
  -16.00000000000003 - 6.6951290300828266e-18im
 -14.000000000000016 + 1.3094008837548432e-15im
 -14.000000000000012 + 2.78841710585368e-17im
 -12.000000000000014 - 2.796888294211598e-17im
 -12.000000000000005 + 7.040218273262716e-16im
 -10.000000000000005 + 2.315044958024994e-16im
  -9.999999999999998 + 1.5382714567093696e-17im
  -8.000000000000004 - 2.1014552699603364e-16im
  -8.000000000000002 + 1.5056230642491578e-17im

Is that what you want? @GiggleLiu

It seems not all minimum eigenvalues are resolved. There should be a 4 fold degeneracy.

giving a smaller tolerance seems to fix that, the default tolerance is sqrt(eps(Float64))

julia> decomp, history = partialschur(-mat(h),nev=10,which=:SR,tol=eps(Float64))
julia> decomp.eigenvalues
10-element Vector{ComplexF64}:
 -16.000000000000057 + 9.443067557412728e-16im
 -16.000000000000025 - 8.669250062378378e-16im
 -15.999999999999991 + 9.371613795987999e-17im
 -15.999999999999973 + 2.1420403513366442e-16im
 -14.000000000000103 + 8.87374359267708e-16im
 -14.000000000000096 - 3.807368843242745e-15im
 -14.000000000000094 - 1.8643396423708125e-15im
 -14.000000000000009 - 1.6927771517392586e-16im
  -13.99999999999999 - 1.6784039256327937e-15im
 -12.000000000000073 + 7.913181027740124e-16im

[hz-xiaxz]

Oops, degeneracy around -14 seems weird...

[hz-xiaxz]

Krylovkit.jl also works well at small tolerance. However, we see instability in the first excited states. I agree block lanczos is neccessary.

julia> using KrylovKit
julia> eigsolve(-mat(h), 10, :SR,tol=eps(Float64))[1]
12-element Vector{Float64}:
 -16.000000000000014
 -16.000000000000007
 -15.99999999999999
 -15.999999999999979
 -14.000000000000057
 -14.000000000000027
 -14.000000000000023
 -14.000000000000014
 -14.000000000000007
 -14.000000000000007
 -13.999999999999988
 -13.999999999999961