We are trying to generate a MPO from OpSum. However, we notice that this process takes much longer time in Julia than in C++. The operator that we are looking at corresponds to a long range Hamiltonian in a spin-1/2 fermion chain. For the same operator string, the MPO(os, sites) function takes more than 11 hours in Julia but only 36 seconds in C++. We are wondering if there is anything we can do to speed up the autoMPO process in Julia.

Thanks for the report. The algorithms are the same in the Julia and C++ versions, so likely we are missing some basic optimizations in the Julia version. Could you share a code that reproduces those timings so we can track this issue?

As a practical workaround, you could try saving the MPO constructed from C++ and load it into the Julia version, if you are interested in using the Julia DMRG/MPS code with those Hamiltonians.

Also @corbett5 has been working on a faster OpSum to MPO conversion function in the Julia version, see GitHub - ITensor/ITensorMPOConstruction.jl: Alternative backends for constructing MPOs from sums of operators., you could try his code.

@mankai_chow even if you arenâ€™t interested in trying out my new package I would be very interested in seeing how the new algorithm performs on your MPO. If youâ€™re willing to share your Julia code (or even just the description of the Hamiltonian) Iâ€™d be happy to give it a spin and let you know how it goes.

Dear Fishman and Corbett :

Thanks a lot for the help.

I cannot upload the file here, so please find a minimal example at the GitHub link https://github.com/mankai-chow/gen_mpo_min.git. On the same computer, the time for the toMPO function in the C++ code is 3.4s and the time for the MPO function in the Julia code is 65s. Since there is precompile issue in Julia, I run the function twice in Julia and take the latter time.

The MPO that we are interested in contains long range interaction. The example has maximal link dimension ~500, but in practice for large system sizes the link dimension could be very large. We overload the Electron type so that we can add an extra quantum number.

I have some trouble getting the new package to work for our model. The new function does not allow that a site has more than one operator acting on it in a term, while our model may need to act terms such as â€śCdagupâ€ť, â€śCdagdnâ€ť simultaneously on the same site. I am not sure how to implement that.

Best regards,

Zheng Zhou ĺ‘¨ćŁ.

Thanks Zheng,

We will be very interested to benchmark your case (Ben, thanks for taking that on). Iâ€™m certain we can get the speed in the Julia version to be the same or even better than the C++ version. This will be very helpful for testing that claim.

Regarding the `"Cdagup", "Cdagdn"`

question, you should be able to put them on the same site like this:

```
os .+= "Cdagup",i,"Cdagdn",i
```

where `i`

is the site and of course you can put a coefficient in front too which I have just left out. If you try that and it leads to an error or incorrect result, please let us know.

Dear Miles :

Thank you for helping. I have checked that this indeed produce an error. My minimal example is

```
os = OpSum()
os += "Cdagup", 2, "Cdagdn", 2, "Cup", 3, "Cdn", 1
hmt = MPO_new(os, sites ; cutoff = 1E-13)
```

and the error message is

```
ERROR: LoadError: A site has more than one operator acting on it in a term.
```

Best regards,

Zheng Zhou ĺ‘¨ćŁ.

This is a self imposed limitation (to achieve greater performance for the case I cared about) of my MPO construction algorithm. However if you provide an extra key-word argument to `MPO_new`

it will convert the `OpSum`

above into `"Cdagup * Cdagdn", 2, "Cup", 3, "Cdn", 1`

which is perfectly valid.

@mankai_chow I created a PR here Changes to work with ITensorMPOConstruction by corbett5 Â· Pull Request #1 Â· mankai-chow/gen_mpo_min Â· GitHub

After my modest changes I can construct the MPO in 0.65 seconds on my laptop.

Dear Corbett :

I have tried the version you provide, and now the AutoMPO process takes 0.73s for me on the same computer, even faster than the C++ version of the code ! That exactly solves our questions. Thank you very much for your help.

Best regards,

Zheng Zhou ĺ‘¨ćŁ.

@mankai_chow what value of `nm`

was the C++ code able to construct in 36 seconds? I am curious because I havenâ€™t compared my algorithm with the ITensorâ€™s C++ implementation before.

I was able to construct a MPO with `nm = 100 and `

tol=4e-8` resulting in an MPO of bound dimension 1788 with the following timings:

```
calc_int_ps: 6.427996 seconds (119.91 M allocations: 3.640 GiB, 18.18% gc time)
constructing opsum: 1543.176156 seconds (165.25 M allocations: 9.714 TiB, 15.88% gc time)
constructing MPO: 51.784839 seconds (71.73 M allocations: 249.206 GiB, 30.54% gc time)
```

After I switched to construting a `ITensorMPOConstruction.OpIDSum`

instead of an `OpSum`

I got the following timings:

```
calc_int_ps: 6.213813 seconds (119.91 M allocations: 3.640 GiB, 18.60% gc time)
Constructing opsum: 6.255719 seconds (119.91 M allocations: 3.645 GiB, 18.68% gc time)
constructing MPO: 46.975967 seconds (55.70 M allocations: 246.379 GiB, 32.90% gc time)
```

I pushed the changes to the PR. You are the first person other than me to use `ITensorsMPOConstruction`

. I am very glad itâ€™s working for you. Let me know if you have more questions or problems.

Dear Corbett :

The case I was mentionning was for n_m=48 (I think this should be the maximal size that DMRG can deal with). We tested this on the Julia version using `OpSum`

and `MPO_new`

and the C++ version. For the AutoMPO process, the Julia version takes 9.3s excluding compilation, and the C++ version takes 39.1s.

Also thank you for the further improved code using `ITensorMPOConstruction`

. Do we need to install any other packages in order to use that ?

Best regards,

Zheng Zhou ĺ‘¨ćŁ.

Thank you for those timings.

OpIDSum is part of ITensorMPOConstruction so you donâ€™t need anything else. In fact the first step of the algorithm is to convert the OpSum to an OpIDSum. But if youâ€™re only going up to 48 sites the OpSum implementation should be plenty fast.

Hi Corbett.

Could you help me have a look at the following error message when using ITensorMPOConstruction ?

```
ERROR: LoadError: MethodError: no method matching set_scalar!(::OpIDSum{ComplexF64}, ::Int64, ::Float64)
Closest candidates are:
set_scalar!(::OpIDSum{C}, ::Integer, ::C) where C
@ ITensorMPOConstruction ~/.julia/packages/ITensorMPOConstruction/VZH0Q/src/OpIDSum.jl:84
Stacktrace:
[1] macro expansion
@ ~/.julia/packages/ITensorMPOConstruction/VZH0Q/src/OpIDSum.jl:158 [inlined]
[2] merge_terms!(os::OpIDSum{ComplexF64})
@ ITensorMPOConstruction ~/.julia/packages/TimerOutputs/Lw5SP/src/TimerOutput.jl:237
[3] macro expansion
@ ~/.julia/packages/ITensorMPOConstruction/VZH0Q/src/OpIDSum.jl:309 [inlined]
[4] prepare_opID_sum!(os::OpIDSum{ComplexF64}, sites::Vector{Index{Int64}}, opCacheVec::Vector{Vector{OpInfo}}, basisOpCacheVec::Vector{Vector{OpInfo}})
@ ITensorMPOConstruction ~/.julia/packages/TimerOutputs/Lw5SP/src/TimerOutput.jl:237
[5] MPO_new(ValType::Type{Float64}, os::OpIDSum{ComplexF64}, sites::Vector{Index{Int64}}, opCacheVec::Vector{Vector{OpInfo}}; basisOpCacheVec::Vector{Vector{OpInfo}}, tol::Int64)
@ ITensorMPOConstruction ~/.julia/packages/ITensorMPOConstruction/VZH0Q/src/MPOConstruction.jl:255
[6] MPO_new
@ ~/.julia/packages/ITensorMPOConstruction/VZH0Q/src/MPOConstruction.jl:245 [inlined]
[7] MPO_new(os::OpIDSum{ComplexF64}, sites::Vector{Index{Int64}}, opCacheVec::Vector{Vector{OpInfo}}; basisOpCacheVec::Vector{Vector{OpInfo}}, tol::Int64)
@ ITensorMPOConstruction ~/.julia/packages/ITensorMPOConstruction/VZH0Q/src/MPOConstruction.jl:264
[8] MPO_new(os::Sum{Scaled{ComplexF64, Prod{Op}}}, sites::Vector{Index{Int64}}; basisOpCacheVec::Vector{Vector{OpInfo}}, tol::Int64)
@ ITensorMPOConstruction ~/.julia/packages/ITensorMPOConstruction/VZH0Q/src/MPOConstruction.jl:284
[9] top-level scope
```

The minimal example I have is

```
using ITensors
using ITensorMPOConstruction
os = OpSum()
os += "Cdag", 2, "C", 2, "Cdag", 1, "C", 1
os += "Cdag", 1, "C", 1, "Cdag", 2, "C", 2
sites = siteinds("Fermion", 2)
operatorNames = [ "I", "C", "Cdag", "N" ]
opCacheVec = [ [OpInfo(ITensors.Op(name, n), sites[n]) for name in operatorNames] for n in eachindex(sites) ]
mpo = MPO_new(os, sites ; basisOpCacheVec = opCacheVec)
```

The error would not occur if either of the two terms were removed. Thank you in advance for your help.

Best regards,

Zheng Zhou ĺ‘¨ćŁ.

Can you start a new post with your new question and/or raise an issue on Github (Issues Â· ITensor/ITensorMPOConstruction.jl Â· GitHub)?