Add CTMRG algorithm for layers of PEPO tensors

[sanderdemeyer]

This builds upon #134, where the CTMRG algorithm now includes InfiniteSquareNetworks without a top and bottom layer, where the top physical leg of the top PEPOTensor is contracted with the bottom physical leg of the bottom PEPOTensor. This is based on the type PEPOTraceSandwich (The equivalent of PEPOSandwich in the previous PR).

This PR can be summarised as

• Define the local sandwich type PEPOTraceSandwich
• Add a constructor for an InfiniteSquareNetwork of PEPOTraceSandwiches based on an InfinitePEPO
• Define the function dagger for an InfinitePEPO
• Add a CTMRG algorithm for a PEPOTraceSandwich
• Generalise the current mpotensor function for PEPSSandwich to both PEPOSandwiches and PEPOTraceSandwiches.
• Add a test where the expectation value for a double layer PEPO is compared to the case where this double layer is now squashed into a single layer
• Add a test where the expectation value for a contraction of a PEPO on a PEPS is compared to the case where this application is squashed into a single layer PEPS

Any suggestions on the implementation and naming of things in this PR is very much appreciated. What anyway needs to be done, is

• Add documentation for e.g. dagger
• Settle on a naming convention
• Maybe add an additional test on the different mpotensor functions that now exist

[codecov]

Codecov Report

Attention: Patch coverage is 78.70370% with 46 lines in your changes missing coverage. Please review.

Files with missing lines	Patch %	Lines
src/algorithms/contractions/ctmrg_contractions.jl	72.41%	40 Missing ⚠️
src/networks/local_sandwich.jl	89.65%	6 Missing ⚠️

Files with missing lines	Coverage Δ
src/algorithms/toolbox.jl	98.14% <100.00%> (+0.07%)	⬆️
src/operators/infinitepepo.jl	48.00% <100.00%> (+5.14%)	⬆️
src/networks/local_sandwich.jl	84.26% <89.65%> (+7.95%)	⬆️
src/algorithms/contractions/ctmrg_contractions.jl	61.57% <72.41%> (+3.48%)	⬆️

... 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.

[sanderdemeyer]

I think everything should be fine now. And if @lkdvos is okay with the naming convention, I think we can merge this.

[sanderdemeyer]

After some minor fixes, I think it should now be okay if everyone is fine with the naming conventions and the current tests @lkdvos and @leburgel.

[lkdvos]

If this is the actual norm, why can't you call norm explicitly? Why do we have to manually divide by the norm, is this not all equivalent to expectation_value?

[sanderdemeyer]

I don't think that there is currently a norm defined for a generic InfiniteSquareNetwork apart from the trivial norm of the tensor itself and the specific case of a PEPS. This would also not be well defined in the generic case, which is why I would be okay with keeping the current version. The expectation_value of an operator (either a PEPO or LocalOperator) of a PEPOTraceSandwich could be defined, but this is another PR (see #117).

Writing this down as an expectation_value of an InfinitePartitionFunction would require us to fuse the layers, which would defeat the purpose of the test.

[lkdvos]

I think my main question is why does the private _contract_site function show up here, if the user-facing function is contract_local_tensor? Also, is that not just value of the network? Basically, if using that private function is the only way to get a sensible result for these contractions, there should be an interface for accessing that

[lkdvos]

To answer the other question: yes I would prefer to write this in terms of contract_local_tensor, so that if we ever decide to change the internal private function, we don't also have to update the tests.

[sanderdemeyer]

I've taken care of the new comments. As a summary, I've changed the docstring to the dot notation, reduced the amount of code in the global scope, and changed the tests to compare the network_values of the mpos instead of the tensors themselves.

[lkdvos]

I think my main question is why does the private _contract_site function show up here, if the user-facing function is contract_local_tensor? Also, is that not just value of the network? Basically, if using that private function is the only way to get a sensible result for these contractions, there should be an interface for accessing that

[lkdvos]

Am I correct that this is making the operator hermitian? If so, I'm not sure the tests for the adjoint are really doing all that much, since the pepo and its adjoint should be equal.

[sanderdemeyer]

I would indeed have liked to omit this line. The problem is that a random PEPO is quite hard to converge. For the (physical) operators that I am computing in the finite-temperature case, they all converge nicely. Since these are not in PEPSKit, I can unfortunately not use them for these tests. For the bosonic case, the convergence is also easier.

I think the tests are still checking the definition of the adjoint, even though it's a hermitian operator (especially in the fermionic case, there are some nontrivial twists). If the _dag function were incorrect, this would just correspond to a different initialisation of the tensor O, and the tests would fail)

[lkdvos]

Is there a major reason you cannot add them to PEPSKit?

[sanderdemeyer]

There reason is that the code is present in my code for cluster expansions, which would be a relatively big addition to PEPSKit (and to be fully correct, requires some functionality for BigFloats that either requires my personal fork of TensorKit (since this is not yet merged, and maybe never will) or, in the long run, maybe MatrixAlgebraKit).

I am still in doubt whether the long-term plan is to add this to PEPSKit, or to keep it as a different package to put on the QuantumKitHub. It will anyway still require a large effort to make everything consistent with PEPSKi.

[lkdvos]

Why not generate the PEPO there and hardcode the entries in the tests here then?

[sanderdemeyer]

That is also a possibility, but this is at least a 2x2x5x5x5x5 fermionic tensor (and I checked, with 184 nonzero elements) and would be quite a hassle to implement. Not to mention that it would be extremely ugly, since the elements are derived from exponentials and are not 'pretty' numbers.

[lkdvos]

I don't really agree with these statements. The hassle you mention is solved with ctrl+c, ctrl+v, and while I agree that there are nicer things, simply copy pasting that into a separate file you include seems very reasonable.

It is just that the longer I look at this PR and the implementations, the more I seem to have the feeling that there are a lot of very brittle things and things that seem to be having wrong twists, and I am not convinced the tests are capturing any of that.

Some alternative suggestions that might be reasonable as well:

• you can map between an infinite PEPS and a PEPOTrace by "bending" the legs downwards and fusing them. This way, we would have more things to cross-check
• at least one physical result should be reproducible in the tests. An energy of a finite temperature state, or something similar, even for a free fermion model would already be a much more convincing statement than I have a random peps, and after hermitianizing and taking its dagger the network value doesn't change.

[lkdvos]

I think here again you have the mathematicians convention, would it not make more sense if you measure O instead of O'? Also, since you explicitly made O hermitian I'm not sure if this is a very good test

[sanderdemeyer]

I don't really get what the problem is with the mathematical convention instead of the physical one. Isn't this just a convention change from $$O$$ to $$O^\dagger$$? Again the same comment about the convergence for $$O^\dagger$$ holds...

[lkdvos]

As the name says, it is definitely a convention, but I would say that it makes a lot more sense to use the physics convention given that the application is physics, or is that crazy? I'm happy to hear if there is a good reason to prefer it the other way around, but it does seem nice to remain consistent.

[sanderdemeyer]

Well, I definitely prefer the mathematical convention, since that is how I learned about adjoints, and so the 'physics convention' seems a bit weird to me. I was also not aware that there were even different conventions, because it seems to me the convention only applies to how you would write the defining equation for the adjoint, and has no influence on the definition itself.

At first glance, I also didn't find anything in PEPSKit or TensorKit where we already chose a convention for this. Is this something of Julia itself?

[lkdvos]

As discussed before, yes, this is the convention of Julia, see dot.

[lkdvos]

If I'm understanding this correct, the networks here are: dot(psi, O * O', psi), dot(M * psi, O * psi), dot(O * psi, O * psi), but I'm somewhat surprised that the contractions below line up at all. Why is M * psi the same as the contract_local_tensor? Shouldn't M * psi be replacing a full layer of PEPO tensors instead, or am I missing things?

[sanderdemeyer]

The M here really is just a local tensor, not a PEPO. In this test, I contract the sandwich of the (infinite layers of) $$\psi$$, $$O$$, $$O^\dagger$$ and $$\overline{\psi}$$, in which I replace $$O$$ with $$M$$ at one site. This is then divided by the network without this replacement. I added this test to make sure that environment of the sandwich is correctly implemented, even if there are multiple PEPO layers and to make sure the code can deal with insertion of operators at specific places (x,y,z coordinates). It is also again a check on whether the adjoint of the PEPO is correct.

[lkdvos]

But why is InfinitePEPS(Mpsi) the same as replacing O with M at one site? That would amount to placing an M at every site?

[sanderdemeyer]

It is not, InfinitePEPS(Mpsi) would indeed correspond to replacing M at every site. That is why this is not used in the calculation of an environment. We only need this because the function network_value requires a InfinitePEPS or InfiniteSquareNetwork and not just a single tensor.

[lkdvos]

Oh I see now. I'm not sure if this is a great way to implement this though, since this very much depends on the exact way that the internal machinery of network_value works. If there is no other way, perhaps it's also okay to at least at a very big warning/disclaimer in the comments, because this is very non-obvious that you are making use of that hack.

The main problem is that in large codebases with multiple people, you really have to be more consistent with what you are expecting from the functions, and network_value is only guaranteed to return you the network value, so if someone else, including yourself in the future forgets that this implemetation detail is being depended on, we would create issues.

[sanderdemeyer]

About your comment on the accessibility of _contract_site and contract_local_tensor: I completely agree, and that is why I would like to add expectation values of a LocalOperator or a PEPOTensor of a PEPOTraceSandwich as a more user-friendly interface. I would like to do this either after or in PR #117.

For now, I can either keep it like this, or change these lines to the equivalent calculation
`
nrm = PEPSKit.contract_local_tensor((1, 1, 1), O, network, env)

nrm_fused = PEPSKit.contract_local_tensor((1, 1, 1), O2, network_fused, env_fused)
`
I don't really have a preference in this regard.

[lkdvos]

I'm very confused why this would be twisting the hth leg of the fuser. Isn't there at least a +1 missing here, since the first leg would be the fused leg?

Given this, I'm also a bit scared about the tests and their coverage, because that seems to imply this isn't showing up at all?
This is definitely possible, since before I wrote the fusers for a single arrow configuration, but I also had stronger assertions on them, which have been loosened here.

[lkdvos]

I don't really agree with these statements. The hassle you mention is solved with ctrl+c, ctrl+v, and while I agree that there are nicer things, simply copy pasting that into a separate file you include seems very reasonable.

It is just that the longer I look at this PR and the implementations, the more I seem to have the feeling that there are a lot of very brittle things and things that seem to be having wrong twists, and I am not convinced the tests are capturing any of that.

Some alternative suggestions that might be reasonable as well:

• you can map between an infinite PEPS and a PEPOTrace by "bending" the legs downwards and fusing them. This way, we would have more things to cross-check
• at least one physical result should be reproducible in the tests. An energy of a finite temperature state, or something similar, even for a free fermion model would already be a much more convincing statement than I have a random peps, and after hermitianizing and taking its dagger the network value doesn't change.

[lkdvos]

To answer the other question: yes I would prefer to write this in terms of contract_local_tensor, so that if we ever decide to change the internal private function, we don't also have to update the tests.

[lkdvos]

As discussed before, yes, this is the convention of Julia, see dot.

[lkdvos]

Oh I see now. I'm not sure if this is a great way to implement this though, since this very much depends on the exact way that the internal machinery of network_value works. If there is no other way, perhaps it's also okay to at least at a very big warning/disclaimer in the comments, because this is very non-obvious that you are making use of that hack.

The main problem is that in large codebases with multiple people, you really have to be more consistent with what you are expecting from the functions, and network_value is only guaranteed to return you the network value, so if someone else, including yourself in the future forgets that this implemetation detail is being depended on, we would create issues.

[lkdvos]

@ogauthe do you maybe have some reference for some high temperature expansion results we could use to verify this implementation on? Alternatively, do you have any inspiration for any testcases that we could use?

[ogauthe]

Coefficients for high temperatures series expansion for the J1-J2 model can be found in from https://doi.org/10.1103/PhysRevB.67.014416
they also work well for J2=0.

I detailed a bit more in #215

[sanderdemeyer]

To get a fermionic test which is reproducable with the current version of PEPSKit, we could look at the spinless fermion model. As mentioned above, I was first planning to add such a test using code on cluster expansions, but since this is a different package and would then have to be included in separate txt files, this is also not ideal. Therefore, I would opt for using Simple Update to produce some finite-temperature state, and then compare with figure 8a or 10a of this paper.

Another thing that would make this chain of PRs in my opinion a bit more logical, is to first think about how we would implement PR #117. While this can certainly be done in parallel, that PR implements the expectation value of (excited) partition function tensors, which is in principle what an expectation value of a LocalOperator of a density operator (represented as a PEPO) would reduce to. In that way, we could test the physical (fermionic) examples already there without any reference to adjoints or PEPOTraces and the CTMRG algorithm on them. We would only have to define a function trace_out to get an InfinitePartitionFunction from an InfinitePEPO and then perform the existing CTMRG algorithm on that. When we have figured out the other PR, we could revisit this PR and see whether it reproduces the same results when we use the CTMRG algorithms on PEPOTraces directly.

On my fork of PEPSKit, I have had a first go at implementing the things above. In this file, it shows how the expectation values reduce to the evaluation of a network of (excited) partition function tensors like in PR #117, while in this and this file the implementation can be found (via the functions expectation_value and contract_local_tensors). The documentation is certainly not yet there, and we could definitely discuss the preferred syntax and implementation, but I think this is a strategy that overcomes the limitation of not being able to use symmetric tensors in PR #117 , and can serve as a building block on which we can solve the problems with this PR.

[lkdvos]

Therefore, I would opt for using Simple Update to produce some finite-temperature state, and then compare with figure 8a or 10a of this paper.

Sounds like a solid plan. I read through that paper, it seems like they are just using the fact that it is free fermions to obtain exact expressions and compare to that. This sounds like an ideal case to put in our tests, since it isn't trivial but also we are certain of the desired output.

For the expectation values, while I get your point, it is a bit strange to tackle that first since expectation_value(::InfinitePEPO, op, ::CTMRGEnv) really needs a CTMRGEnv, which we cannot obtain without first implementing the CTMRG.
Am I understanding it correctly that for your implementation, since you are converting to a partition function, you would have to provide the CTMRGEnv for the partition function, which only coincides when there is a single layer?

I think that actually the reduced density matrix approach might be a better fit after all. It should be somewhat straightforward to generate the contractions for building these reduced density matrices, after which applying operators really is just tr(O * rho) / tr(rho). This turns out to be a lot cheaper as well, since we don't have to perform two contractions (numerator and denominator), and I verified that the leading cost of computing the density matrix is actually equal to the leading cost of immediately applying an operator, so there really is no lost performance.

I do like the idea as well of having additional consistency checks between the different algorithms: For example turning the density matrix into a PEPS by fusing the legs, we can evaluate expectation values that correspond to the partition function at 2beta. Similarly, the network value (partition function) should be the same for the partition function, the PEPS, and the PEPO approach, which provides good checks.

I'm not entirely convinced that the expectation_value of a partition function is something we want to throw in the mix already, mostly because it seems like the interface still is missing something. It would be great if we could formulate it in a similar framework (density-matrix-like, so we would obtain an object we can then use to evaluate arbitrary operators), but I would have to think through how (and if?) that translates well. It sounds a lot like source terms etc, but I don't have that knowledge directly available.

[sanderdemeyer]

The point I wanted to make was that the CTMRG algorithm for partition functions is already implemented. You could indeed write is such that you need a CTMRGEnv of the partition function. Maybe a more natural option is to require a envspace, so we can call something like env0 = CTMRGEnv(network_pf, envspace). Then you only need an environment bond dimension. The approach indeed only works for 1 layer, which would be fine for me. I don't think many people want to represent a density operator as a multi-layer PEPO. If we do want to allow for multiple layers, we could still convert to an InfinitePartitionFunction by fusing the legs.

What do you mean exactly with 'the density-matrix approach'? Do you mean the implementation of this PR? Like you said, this reduces to tr(O * rho) / tr(rho), which is indeed just a contraction of an InfinitePartitionFunction (which was the point I was trying to make). Or am I misunderstanding you?

The current interface of the expectation_value of an InfiniteParititionFunction is indeed not optimal. I would just like to add it since expectation values of PEPOs naturally reduce to them, and since this framework would also allow us to generalize this to n-site operators in a rather straightforward manner. Whether we should export it and whether we should call it an expectation_value of an InfinitePartitionFunction is certainly up for discussion.

[lkdvos]

What do you mean exactly with 'the density-matrix approach'?

See #230 : instead of inserting the operators immediately, we just literally compute the reduced density matrix and then apply the operator. This is generalizable to PEPOtrace, and retains the same interface.

Maybe a more natural option is to require a envspace, so we can call something like env0 = CTMRGEnv(network_pf, envspace). Then you only need an environment bond dimension.

I think my main point was that you start with expectation_value(::InfinitePEPO, O, env), and I can't see what env is here. Simply providing spaces is not enough, you'd still have to converge it through some algorithm, so you would have to already converge the environment for the partition function, which then defeats the purpose a bit.

The approach indeed only works for 1 layer, which would be fine for me. I don't think many people want to represent a density operator as a multi-layer PEPO.

If the approach only works for a single layer, then really we should just use partition functions to begin with, since it is more efficient to trace the pepo once instead of keeping the legs separate. This is a viable strategy btw, if that makes things easier. Similarly, for CTMRG and a double layer PEPO, it is actually easier to map this to an InfinitePEPS double layer, since the contractions would be exactly the same but with slightly less overhead because of a reduced number of legs. (This is not true for the reduced density matrix of a double layer, since there we really want to trace the legs as soon as possible).

Also note that in the context of purification with a hermitian and (possibly real, or time-inversion invariant, I forgot) operator, it is actually beneficial to write the reduced density matrix in terms of a double layer object since that would become positive definite by construction, which might give you more physical results: rho(beta) = tr(rho(beta)/2 * adjoint(rho(beta/2)).

The current interface of the expectation_value of an InfiniteParititionFunction is indeed not optimal. I would just like to add it since expectation values of PEPOs naturally reduce to them

This is not entirely true (anymore?) since it's actually more convenient to work with reduced density matrices (this was pointed out by @ogauthe in the linked PR), which works well for PEPOs but not for partition functions.