[Discussion] Finite temperature imaginary time evolution

I plan to add the feature to calculate ρ = exp(-βH) via simple/full update. The idea is that we start from a product state with an "ancilla" space attached to each PEPS tensor:

$$|\psi(0)\rangle = \prod_j \sum_s {|s^{phy}_j, s^{anc}_j\rangle}$$

Then we perform imaginary time evolution on the physical space only (i.e. we apply identity gates to the ancilla space)

$$|\psi(\beta)\rangle = \exp(-\frac{1}{2}\beta H) {|\psi(0)\rangle}$$

Finally, exp(-βH) is obtained by tracing out the ancilla

$$\rho(\beta) = \exp(-\beta H) = \mathsf{tr}_{anc} {|\psi(\beta)\rangle} {\langle\psi(\beta)|}$$

The issue is whether it is worthwhile to keep the physical and the ancilla spaces separated during the evolution, which will need a new struct similar to InfinitePEPS but with more than one codomain spaces, and duplicate lots of code.

A compromised version is we only split the physical and the ancilla spaces when applying the gate. Then they are fused for the rest of an SU/FU step. But from what I know, TensorKit does not provide a "efficient" way to split/fuse indices of a tensor (without contraction with a fuser tensor). So maybe it won't affect performance too much if we just always keep the physical and the ancilla spaces fused all the time (this will require the ancilla space to be fused into the gates, LocalOperators etc)?

@ogauthe

I took a look at the TeNeS package and the accompanied paper. There they briefly discussed two approaches:

Evolving the purified state ${|\psi\rangle}$ as explained above. Theoretically the PEPO of ρ(β) can be obtained by tracing out the ancilla leg, but the resulting bond dimension will become D^2, where D is the bond dimension of ${|\psi\rangle}$, which is non-practical. But the advantage is that ρ(β) is guaranteed to be positive definite. PhysRevB.99.035115 further propose that one can apply a "disentangler" (a unitary operator on the ancilla spaces) to "reduce the required D".
Directly evolving ρ(β) = exp(-βH) as a PEPO, with the Trotter gate applying on both sides ρ(β) = exp(-βH / 2) ρ(0) exp(-βH / 2). This allows pushing to a large bond dimension D of the PEPO (though I doubt whether this really has advantage over the D^2 in the purified approach), but the positive definiteness of ρ(β) is no longer guaranteed due to approximations made during the bond truncation. This is the approach used in TeNeS.

@ogauthe Which approach do you recommend to be implemented in PEPSKit?

@sanderdemeyer

@sanderdemeyer I also wonder how SU/FU complement the cluster expansion method you are developing.

One thing I noticed in your and @ogauthe's code is that you choose the ancilla space to be the dual of the physical space. But intuitively, the purification idea should lead to ancilla space = physical space, without the need to take the dual?

I don't think these approaches are actually all that different, and they just constitute different algorithms for the exact same thing. The starting point is mostly just a density matrix $\ket{\psi}\bra{\psi}$, typically the infinite temperature one where this is effectively the identity operator.
We can represent this object as a PEPO (although the infinite temperature one is a "product PEPO"), or a PEPS with two physical legs. These representations are equivalent, and we can go between the two representations at any point, simply by bending the PEPO leg down, or the PEPS leg up. (See also how mapping an MPO to an MPS with two physical legs has that same idea).

Then, we can use imaginary time-evolution to evolve the temperature, typically in one of the following ways:

We apply trotter gates on one of the physical legs (note that this can be done both in the PEPS as well as the PEPO representation)
We apply trotter gates on both of the physical legs
We apply "evolution operators" to decrease the temperature (i.e. $\rho(2\beta) = \rho(\beta) * \rho(\beta)$

Here, clearly (1) manifestly breaks the bra-ket (hermitian) symmetry of the density matrix, and (3) tends to be very expensive since it involves triple layer environments.

Then, in order to evaluate expectation values, we require $tr(O\rho)$. This can again be evaluated in different ways:

Explicitly trace out the ancilla, obtain a partition function and compute CTMRG environments (single layer), then compute the remaining local contraction
Rewrite $tr(O\rho(\beta)) = tr(\rho(\beta/2) O \rho(\beta/2)) = dot(\rho(\beta/2), O, \rho(\beta/2))$, where the last equality holds only for real density matrices. Then contract the double layer object by fusing physical leg and ancilla, doing regular peps CTMRG and computing the remaining local contraction.

Here the first one scales as $D$, the second as $D^2$, but the second only requires $\beta / 2$, and is manifestly hermitian again.

Note also that we can "mix-and-match" between the different methods at any point.
So practically speaking, it doesn't really matter which data structure we choose, it is simply which one we find most convenient and a convention.

All of this also gives an idea about why, in the interpretation of the PEPS with two physical legs, one of them is dual. This stems from the idea that for a purified wavefunction, you get $\rho = \ket{\psi}\bra{\psi}$, and not $\bra{\psi}\ket{\psi}$. Because of the bra-ket, when you then reinterpret both of them as physical spaces for the PEPS (bending the legs), one of them has the opposite arrow as the other. Note also that this has to be the case if you want to be able to trace over these legs, otherwise the arrows simply don't match.
In principle, it's possible to then do a local transformation on one of the legs to flip the arrow, to obtain two arrows in the same direction, but that would still be a space that is isomorphic to the dual space, not the regular one.

I agree with the statements made in the TeNeS paper and therefore have a preference for the density operator w.r.t. the purified state. I do get the appeal for the purification, and certainly if we would implement a method for disentangling, this could maybe be the better choice in certain cases. Therefore, it may be best to support both, since the only differences between the two choices would be the construction of the gates and the evaluation of expectation values. In that regard, fusing the legs into a PEPS in the purified case seems the most convenient to me. For the expectation values, in the case of a density operator, I just trace out the physical leg to calculate the environment of the partition function. In the purified case, I use the last equation of Lukas's 2nd point, i.e. $$<\rho(\beta / 2) | O | \rho(\beta / 2)>$$, which is just the expectation value of an operator (with fusers) of a PEPS.

I have used SU to do the time evolution in the purified case, but have noticed that it is way less accurate than other approaches, even compared to methods that are computationally quite cheap. The current state-of-the-art is really PEPO-PEPO application and truncation. Most people use the Trotter decomposition, but indeed we now use cluster expansions (see e.g. this paper). Another paper showing that this method is more accurate than Trotter and SU is hopefully coming soon. The truncation methods that I use for this are in currently in the cluster expansions package. When they are rigorously tested and people are interested, this could be merged into PEPSKit. Even though I would very much like to know how well FU performs with respect to SU, I doubt that it would be better than the other methods. So right now, I am not planning to do any more finite-temperature simulations using SU.

I was using a non-dualed ancilla space based on the following interpretation:

When the ancilla is a dual space, I think you're actually already evolving the PEPO exp(-βH) itself, instead of its purification. In this case the ancilla is the "physical domain" of the PEPO, instead of an additional physical space of the ancilla object.

I don't think this interpretation requires the ancilla space to be non-dual though?

I think one intuitive argument for using a dual one is that if you want the infinite temperature state to be representable as a product state, you need the physical space and the ancilla to fuse together to the vacuum. (vac in P * A). If you then consider for example a bosonic physical space P = U1Space(i => 1 for i in 0:cutoff), this means that the ancilla, if the arrow is pointing in the direction on the left of your picture, has to now be A = U1Space(-i => 1 for i in 0:cutoff) (which is equal to flip(P), not dual(P)). On the flip side (pun intended), if the arrow points in the other direction, you get A = dual(U1Space(i => 1 for i in 0:cutoff)).

Obviously this is all just a matter of convention, but I think I like this picture since then at least one obvious way of purifying the density matrix is simply "bending the input pepo physical space down", i.e. choosing A = dual(P), in which case the two interpretations match up.
$\rho \colon P \mapsto P \qquad \simeq \qquad \psi \colon \mathbb{C} \mapsto P \otimes P^*$

Of course, in principle none of the code HAS to depend on this convention, but as we have noticed fermions do complicate that slightly.

@lkdvos

@lkdvos Do we already have something like expectation_value(ρ::InfinitePEPO, O::LocalOperator, env::CTMRGEnv) (where env is the environment of the partition function tr(ρ)) to calculate tr(ρO) / tr(ρ)?

I don't think so, although many of the components might already be there. I think the main thing to support is the reduced density matrices, after which the rest of the code should be the same as for regular PEPS

@sanderdemeyer

@sanderdemeyer In cluster expansion approach, is the expectation values also obtained from tr(ρO)/tr(ρ) using CTMRG for partition functions after constructing the PEPO form of ρ? I took a quick look at your package, and found that the list of observables is one of the inputs to the time evolution, and there are many VOMPS references, which confuses me a bit...

That is indeed how I calculate the observables. The alternative is to use a boundary-MPS approach (with VUMPS from MPSKit) to calculate local observables. VOMPS refers to the truncation scheme (in contrast to the files local_truncation and global_truncation, maybe I should rename the vomps file to variational_truncation) after applying a PEPO to another PEPO. This is currently called VOPEPO, since this is the PEPO version of what is implemented in MPSKit for MPS (see this file), which approximates the application of an MPO to an MPS as a new MPS, whereas here we approximate the application of a PEPO on another PEPO as a new PEPO.

@sanderdemeyer

@sanderdemeyer If you have time, can you help make a PR on calculating tr(ρO)/tr(ρ) with CTMRG first?

Currently we have reduced_densitymatrix(inds, ket, bra, env). Then (as already outlined by @lkdvos) if we similarly define reduced_densitymatrix(inds, ρ::InfinitePEPO), then adding MPSKit.expectation_value(ρ::InfinitePEPO, O, env) should require minimal modification from the existing MPSKit.expectation_value(peps, O, env).

In the mean time I can work on the 3-site SU for finite temperature.

I am definitely interested in this, but I'm going on a vacation starting at the end of this week, and I will only be back on the 12th of October. So I can definitely help, but from next week onwards, I can only comment on some decisions that you make, but I can't really help with the development itself.

If I find time tomorrow I'll just quickly add it in, it's just writing a couple contractions that are almost all already there

Just back from vacations. I will have a look tomorrow.

I now have an idea on next-nearest neighbor (NNN) full update. Instead of alternatively optimizing in a 3-site environment (like https://arxiv.org/abs/1711.07584), we can do a further (not or slightly truncated) SVD in the middle of the 3-site to reduce it to two nearest neighbor gates:

        ↓       ↓
        g2′-←--g3
        ↓       ↓
g1 -←-- g2
↓       ↓

The dimension of the physical space of the middle site will temporarily grow. The advantage is that NNN update can now reuse NN code. If this works well also for SU, the OBC-MPS truncation code currently written for 3-site SU can then be removed.

Just trying to wrap my head around this idea, can I think of this as a somewhat generalized/altered swapgate approach? i.e., for a purely next-nearest neighbour interaction (so not a 3-body term), this looks a lot like first doing a swap gate, then applying the NN gate, but instead of having to swap the sites back at the end you factor the gate slightly differently to avoid having to do a final swap gate?

The first thing I wonder is if we would get into trouble that the "physical leg" between the two gates is now altered, since it depends on the gate decomposition.
The second thing is that I guess if we make the code flexible enough, this more or less boils down to being able to do arbitrary Trotterization schemes where all gates are nearest neighbour (and thus might require swap gates).

In any case, looks like a cool idea to pursue, although I'm not sure if removing the other code is fair, as there might be cases where this works better?

a somewhat generalized/altered swapgate approach

In current 3-site SU implementation, the 3-site MPO gate is constructed as follows:

Take tensor product of the gate on sites 1 & 3 and the identity operator on site 2, so a 3-site gate is obtained.
Permute the indices to reorder the sites from (1, 3, 2) to (1, 2, 3). I think you may say swap gates are indeed involved here.
Decompose the 3-site gate to MPO.

if we would get into trouble that the "physical leg" between the two gates is now altered

This is an intermediate step if we only decompose the 3-site MPO into 2-layers of NN gates, which can be regarded as a price we pay to reduce the number of gate layers from 3 (swap gate approach) to 2.

make the code flexible enough

This requires a dedicated struct to handle the Trotter gates. I can't think that far yet...

removing the other code

For SU I may keep the old code and dispatch on the type of gates. But for FU I surely will only implement the new approach.

Navigation Menu

Repository navigation

[Discussion] Finite temperature imaginary time evolution #247

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

[Discussion] Finite temperature imaginary time evolution #247

Description

Activity

Yue-Zhengyuan commented on Aug 27, 2025

Yue-Zhengyuan commented on Aug 27, 2025

lkdvos commented on Aug 27, 2025

sanderdemeyer commented on Aug 27, 2025

Yue-Zhengyuan commented on Aug 27, 2025

lkdvos commented on Aug 27, 2025

Yue-Zhengyuan commented on Aug 28, 2025

lkdvos commented on Aug 28, 2025

Yue-Zhengyuan commented on Sep 1, 2025

sanderdemeyer commented on Sep 1, 2025

Yue-Zhengyuan commented on Sep 1, 2025

sanderdemeyer commented on Sep 1, 2025

lkdvos commented on Sep 1, 2025

ogauthe commented on Sep 3, 2025

Yue-Zhengyuan commented on Nov 11, 2025

lkdvos commented on Nov 11, 2025

Yue-Zhengyuan commented on Nov 11, 2025

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions