Exact MPO option for the julia version

Question

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

Answer 1

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

Comments

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!

---

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

---

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.

---

Thanks. I can reproduce it now & will look into it.

---

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!