+1 vote
asked by (170 points)

Hi,

I have been looking for a way to create an exact MPO when using the julia version of ITensor. I know that in the C++ version, one can use the function toMPO with the argument "Exact", but I could not find a way to do this in the julia version. Is there already a way to do it? And if not, would it be possible to add one?

Best,

Rafael

1 Answer

0 votes
answered by (49k points)

Hi, good question but we don't offer the "Exact" back-end for AutoMPO right now in the Julia version. We just haven't ported it over, since the other backend (which I guess could be called "Compressing" or "SVD") is generally more useful as it can handle a much wider set of inputs.

However, it's important to know that the current backend is very accurate. It does make some approximations but these are basically to truncate singular values of certain matrices going into the MPO which are below 1E-13, so machine precision zero essentially.

Finally, you can test the quality of the MPO that is produced by computing matrix elements with product state MPS. You can use the productMPS function to make these MPS, then the inner(psi,H,phi) function to compute the matrix elements, to check whether they match the exact result you expect for your Hamiltonian.

Please let me know if you have any follow up questions about any of this.

Miles

commented by (170 points)
Hi Miles,

Thanks for your reply. The reason I asked for this feature is that it seems like the MPO approximation seems to sometimes trigger an error saying the Hamiltonian does not appear to be hermitian, when it actually should be.

I wrote some minimal code to exemplify the problem. It creates a Hamiltonian for a spinless fermionic system of 4 sites arranged in a square, with first and second-neighbor hopping. The second neighbor hopping is complex.

    using ITensors
    
    let
    
      N = 4
      sites = siteinds("Fermion",N,conserve_qns=true)
      t1 = 1.0
      t2p = 1.0+0.1*im
      t2n = conj(t2p)
        
      ampo = AutoMPO()
    
      #vertical bonds
      for j in [1,3]
        ampo += -t1,"Cdag",j,"C",j+1
        ampo += -t1,"Cdag",j+1,"C",j
      end
    
      #horizontal bonds
      for j in [1,2]
        ampo += -t1,"Cdag",j,"C",j+2
        ampo += -t1,"Cdag",j+2,"C",j
      end
    
      #second-neighbor bonds
      for j=1
        ampo += -t2p,"Cdag",j,"C",j+3
        ampo += -t2n,"Cdag",j+3,"C",j
      end
    
      for j=2
        ampo += -t2p,"Cdag",j,"C",j+1
        ampo += -t2n,"Cdag",j+1,"C",j
      end
    
      H = MPO(ampo,sites)
      println(ampo)
      sweeps = Sweeps(5)
      maxdim!(sweeps,10,20,100,100,200)
    
      init_state = ["Occ","Emp","Emp","Emp"]
    
      psi0 = productMPS(sites,init_state)
    
      energy,psi = dmrg(H,psi0,sweeps)
    
      return
    end

When I run this code, most of the time it works fine, but sometimes, without changing any lines of code, after printing the MPO it will throw the following error:

LoadError: operator does not appear to be hermitian: NaN vs NaN
Stacktrace:
 [1] error(::String) at ./error.jl:33
 [2] initialize(::KrylovKit.LanczosIterator{ProjMPO,ITensor{3},KrylovKit.ModifiedGramSchmidt2}; verbosity::Int64) at /home/rafael/.julia/packages/KrylovKit/oJ59b/src/krylov/lanczos.jl:78
 [3] eigsolve(::ProjMPO, ::ITensor{3}, ::Int64, ::Symbol, ::KrylovKit.Lanczos{KrylovKit.ModifiedGramSchmidt2,Float64}) at /home/rafael/.julia/packages/KrylovKit/oJ59b/src/eigsolve/lanczos.jl:10
 [4] eigsolve(::ProjMPO, ::ITensor{3}, ::Int64, ::Symbol; kwargs::Base.Iterators.Pairs{Symbol,Real,NTuple{4,Symbol},NamedTuple{(:ishermitian, :tol, :krylovdim, :maxiter),Tuple{Bool,Float64,Int64,Int64}}}) at /home/rafael/.julia/packages/KrylovKit/oJ59b/src/eigsolve/eigsolve.jl:163
 [5] macro expansion at /home/rafael/.julia/packages/ITensors/ka3Bv/src/mps/dmrg.jl:161 [inlined]
 [6] macro expansion at /home/rafael/.julia/packages/TimerOutputs/dVnaw/src/TimerOutput.jl:206 [inlined]
 [7] macro expansion at /home/rafael/.julia/packages/ITensors/ka3Bv/src/mps/dmrg.jl:160 [inlined]
 [8] macro expansion at ./timing.jl:233 [inlined]
 [9] dmrg(::ProjMPO, ::MPS, ::Sweeps; kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at /home/rafael/.julia/packages/ITensors/ka3Bv/src/mps/dmrg.jl:138
 [10] dmrg at /home/rafael/.julia/packages/ITensors/ka3Bv/src/mps/dmrg.jl:90 [inlined]
 [11] #dmrg#605 at /home/rafael/.julia/packages/ITensors/ka3Bv/src/mps/dmrg.jl:22 [inlined]
 [12] dmrg at /home/rafael/.julia/packages/ITensors/ka3Bv/src/mps/dmrg.jl:21 [inlined]
 [13] top-level scope at /home/rafael/Documents/Haldane/ErrorExample.jl:45
 [14] include_string(::Function, ::Module, ::String, ::String) at ./loading.jl:1088
in expression starting at /home/rafael/Documents/Haldane/ErrorExample.jl:3

Since I am not changing anything on the code, I believe the source of this error is probably the approximation made by the MPO. Do you agree? And if so, do you think adding the option for an exact MPO is the best way to stop this from hapenning? Maybe there is some workaround I haven't thought of yet.

Thanks for the help!
commented by (49k points)
Hi, so thanks for the bug report. It's possible that the problem is the approximation made by AutoMPO, but a more likely origin of the bug is not so much the approximate nature of AutoMPO but just that I might have missed a conj somewhere in that code for complex Hamiltonian cases (so would apply to even the "Exact" backend too maybe).

However, running your code above I didn't get any errors, even after running it a number of times. Also I multiplied all of the MPO tensors together, then took the Hermitian conjugate of that single "H" tensor and it was equal to H, so it appears to be correctly Hermitian.

Could you tell me other steps to reproduce? Does it only happen quite rarely? Also have you upgraded to the latest version of ITensors.jl ?

Best,
Miles
commented by (170 points)
Yes, I am using the latest version of ITensor. For me, I would say the error message comes up approximately once every ten runs. There are no other steps to reproduce it, I simply run the code as is, and occasionally it throws that error message.
commented by (49k points)
Thanks. I can reproduce it now & will look into it.
commented ago by (49k points)
Thanks again for your bug report. We believe this issue is fixed by recent versions of ITensors.jl, probably because they are now using an updated version of KrylovKit.jl which supplies our Lanczos implementation. But please let us know if the bug returns, of course!
Welcome to ITensor Support Q&A, where you can ask questions and receive answers from other members of the community.

Formatting Tips:
  • To format code, indent by four spaces
  • To format inline LaTeX, surround it by @@ on both sides
  • To format LaTeX on its own line, surround it by $$ above and below
  • For LaTeX, it may be necessary to backslash-escape underscore characters to obtain proper formatting. So for example writing \sum\_i to represent a sum over i.
If you cannot register due to firewall issues (e.g. you cannot see the capcha box) please email Miles Stoudenmire to ask for an account.

To report ITensor bugs, please use the issue tracker.

Categories

...