Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Feature request: block lanczos for solving degenerate ground states #87

Open
GiggleLiu opened this issue May 28, 2024 · 7 comments
Open

Comments

@GiggleLiu
Copy link

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
Copy link

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

@Jutho
Copy link
Owner

Jutho commented Jul 17, 2024

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

@hz-xiaxz
Copy link

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
Copy link
Author

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
Copy link

hz-xiaxz commented Aug 24, 2024

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
Copy link

Oops, degeneracy around -14 seems weird...

@hz-xiaxz
Copy link

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

4 participants