Simple update for PEPO

[Yue-Zhengyuan]

This PR adds the ability of simple update to evolve a PEPO. The most well-known application is calculating exp(-βH) at finite temperature.

• Overloaded MPSKit.infinite_temperature_density_matrix to creates the trivial PEPO (i.e. ρ = exp(-βH) at infinite temperature) given a Hamiltonian (as a LocalOperator).

• simpleupdate gets a new kwarg gate_bothsides::Bool, controlling how the evolution of PEPO is done:

  • true: (default) exp(-H dt/2) ρ exp(-H dt/2). This is also how TeNeS does it.
  • false: exp(-H dt) ρ (corresponding to the purified approach, without fusing the physical and the ancilla spaces)
    The total evolution time β is always maxiter * alg.dt. In the former (default) approach the time slice for each Trotter gate is actually dt / 2. This kwarg has no effect for PEPS evolution.

• The kwarg bipartite does not apply to finite temperature evolution.

Things to do in follow-up PRs

• Implement the 3-site version to allow Hamiltonians with next-nearest neighbor terms (e.g. J1-J2; done in Simple update for PEPO (3-site version) #262).
• Add reduced_densitymatrix and expectation_value support for calculating tr(ρΟ) / tr(ρ). (Done in Add double-layer-pepo reduced densitymatrix #252 and Add single-layer-pepo reduced densitymatrix #254)

[codecov]

Codecov Report

❌ Patch coverage is 97.61905% with 2 lines in your changes missing coverage. Please review.

Files with missing lines	Patch %	Lines
src/operators/localoperator.jl	50.00%	2 Missing ⚠️

Files with missing lines	Coverage Δ
src/PEPSKit.jl	100.00% <ø> (ø)
src/algorithms/contractions/localoperator.jl	94.98% <ø> (ø)
src/algorithms/time_evolution/evoltools.jl	97.82% <100.00%> (+0.68%)	⬆️
src/algorithms/time_evolution/simpleupdate.jl	100.00% <100.00%> (ø)
src/environments/suweight.jl	84.21% <100.00%> (+1.51%)	⬆️
src/operators/localoperator.jl	87.71% <50.00%> (-1.38%)	⬇️

... and 1 file with indirect coverage changes

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[Yue-Zhengyuan]

Here are benchmark results on the transverse field Ising model (J = 1 and g = 2), using PEPO D = 8 and CTMRG χ = 16, and evolution time step dt = 0.001. The non-purified approach smears out the transition compared to HOTRG, and produce a smaller x-magnetization in the high-T phase, so it is indeed bad as mentioned by @sanderdemeyer in #247. But the smearing out is much reduced after purification.

I still need to compare with TeNeS (which also uses simple update) to see if they also suffer.

[Yue-Zhengyuan]

Using parameters given in 2501.07777, I can reproduce x-magnetization for transverse field Ising model obtained by TeNeS. They match really well:

Note that they use S = σ/2 in the Hamiltonian, and varied dt during the evolution (the τ in their paper is the Trotter gate time slice, i.e. alg.dt / 2). I haven't compared the energy since we currently don't directly allow more than one "impurity tensor" in the partition function, but I'm confident they will also match well. TeNeS also showed that their SU energy and x-magnetization match the QMC result, though they somehow (maybe deliberately) leave out the z-magnetization, which is the order parameter.

[lkdvos]

Looks great, thanks for the work!

I left some comments throughout, but there are some additional things I'm not fully convinced about that we could discuss:

I'm not super convinced by the gate_side keyword arguments.
It's a bit strange to me to use a symbol at the top level, and then suddenly change it to an integer a bit further down the callstack. In this case, it might make more sense to have either the symbols all the way through, or @enum GateSide left right both if you prefer to work with integers to begin with. There should probably also be some checking that the symbol is one of the allowed symbols to begin with, since this sometimes gets the compiler confused about return types because in principle neither of the branches could be chosen, and this is easy to have subtle type stability bugs.

Additionally, some of these functions to me feel a bit like they are better left as separate functions, rather than using a keyword to distinguish between them (especially the lower down functions). This might be a bit subjective though, feel free to disagree!

Am I correct in assuming that :domain and :codomain side could also be implemented by correctly permuting the PEPO/gate, so really there is only one of them that is relevant?
Additionally, for the :both option, it might be more efficient to really map it to a PEPS and do the update there, as this would require less overhead from the permutations etc? The interesting part is that this would really apply both at the same time, rather than sequentially.

Overall, this definitely looks great and I think the results are as expected from simple update results for finite temperature. It woudl definitely be nice to have the expectation values/reduced density matrices working in the near future as well, and have some tests for fermionic systems, but I'm okay with adding that in later.

[lkdvos]

Does it make sense to have the default gate_side be :both? I might be mistaken but I think that one is the least common approach.

[Yue-Zhengyuan]

In the next fix I shall temporarily set it to :codomain before we finally decide on the interface.

[Yue-Zhengyuan]

Now I changed it to gate_bothsides::Bool, which defaults to true (following TeNeS). A side effect is that :domain mode is removed. For InfinitePEPS I just further set it to false inside the SU function, which does not look very elegant. Is there a better way to make this kwarg applicable to InfinitePEPO only?

[lkdvos]

Note that this function is, if I'm not mistaken, equivalent to

PEPSKit.jl/src/operators/localoperator.jl

Line 192 in df7a2c5

function _fuse_ids(op::AbstractTensorMap{T, S, N, N}, Ps::NTuple{N, S}) where {T, S, N}

[Yue-Zhengyuan]

I added a version of it without Ps. In this case Ps is set to the domain of op (converted to a tuple). I removed attach_ancilla for LocalOperators, though.

[Yue-Zhengyuan]

Am I correct in assuming that :domain and :codomain side could also be implemented by correctly permuting the PEPO/gate, so really there is only one of them that is relevant?

Their difference lies in which physical axis of the PEPO tensor is transferred to the "reduced tensor" when performing QR/LQ in _qr_bond (currently specified by gate_ax there). So in principle :domain is the same as :codomain with a PEPO from dual spaces to dual spaces (which is something I don't feel comfortable with). What do you think if I just remove the :domain mode?

Additionally, for the :both option, it might be more efficient to really map it to a PEPS and do the update there, as this would require less overhead from the permutations etc? The interesting part is that this would really apply both at the same time, rather than sequentially.

I did consider applying the gate on both sides at the same time. But in this way to bond dimension truncation will be from $d^2 D$ (instead of $d D$) back to $D$, which use more memory and may result in a larger truncation error (though I'm not sure about the latter point). Anyway the sequential approach is just simpler to implement, and tensor permutation is never the bottleneck of the algorithm.

[Yue-Zhengyuan]

The test on transverse field Ising model is now updated to use the new expectation_value for PEPO. I think it remains to decide on the interface (in particular, gate_bothsides and its default value) before merging.

[lkdvos]

I think I'm okay with your proposed interface, I can't really come up with something better myself right now either, so I think this is just what it has to be

[lkdvos]

Can you check if the following works in the type domain?