Add CTMRG support for regular 2D partition functions

[leburgel]

Rebased from @Sander-De-Meyer's initial work here.

Some todos left:

• Agree on struct names
• Change indexing convention for local partition function tensor to W ⊗ S ← N ⊗ E and document
• Reorganize partition function code a bit (i.e. try not to mix in with PEPS code in the same files if possible)

[codecov]

Codecov Report

Attention: Patch coverage is 53.70370% with 150 lines in your changes missing coverage. Please review.

Files with missing lines	Patch %	Lines
src/algorithms/contractions/ctmrg_contractions.jl	50.00%	49 Missing ⚠️
src/states/infinitepartitionfunction.jl	24.00%	38 Missing ⚠️
src/networks/infinitesquarenetwork.jl	61.29%	24 Missing ⚠️
src/algorithms/ctmrg/sparse_environments.jl	33.33%	20 Missing ⚠️
src/environments/ctmrg_environments.jl	56.66%	13 Missing ⚠️
src/algorithms/contractions/localoperator.jl	73.33%	4 Missing ⚠️
src/networks/tensors.jl	66.66%	2 Missing ⚠️

Files with missing lines	Coverage Δ
src/PEPSKit.jl	87.50% <ø> (ø)
src/algorithms/ctmrg/ctmrg.jl	87.27% <100.00%> (+0.73%)	⬆️
src/algorithms/ctmrg/sequential.jl	98.27% <ø> (ø)
src/algorithms/toolbox.jl	97.43% <100.00%> (+0.71%)	⬆️
src/operators/infinitepepo.jl	23.80% <100.00%> (+5.62%)	⬆️
src/states/infinitepeps.jl	69.38% <100.00%> (+3.45%)	⬆️
src/states/infiniteweightpeps.jl	93.33% <ø> (ø)
src/networks/tensors.jl	66.66% <66.66%> (ø)
src/algorithms/contractions/localoperator.jl	95.19% <73.33%> (+0.24%)	⬆️
src/environments/ctmrg_environments.jl	61.70% <56.66%> (-1.02%)	⬇️
... and 4 more

[lkdvos]

Starting to look good!
I left some remarks alreayd, but some general remarks here:

For the docstrings, I don't think, it makes sense to have separate docstrings for both the partition function and PEPS (and eventually PEPO) implementations. Rather, we could use A in the figure, and then specify that A might be either bra-ket, partition function or bra-pepo-ket combinations. As a result, we should probably also agree on a style of drawings (let's do this in a separate PR) that uses the same character for drawing tensor connections, and maybe either never draw boxes around tensors or always draw boxes.
(I am assuming that you left some of the docstrings unattached to functions on purposes as this is still a WIP port?)

For naming things, InfinitePartitionFunction seems like the best we can do, since I really cannot come up with anything better. Given that it is rather long, how would you feel about InfinitePF, similar to InfinitePEPS? (Along with PFTensor etc)

[leburgel]

For the docstrings, I don't think, it makes sense to have separate docstrings for both the partition function and PEPS (and eventually PEPO) implementations. Rather, we could use A in the figure, and then specify that A might be either bra-ket, partition function or bra-pepo-ket combinations.

Sounds good, I'll try to standardize as much as I can here. I agree we should be more systematic about the diagrams, I'm happy to put some effort into that soon.

For naming things, InfinitePartitionFunction seems like the best we can do, since I really cannot come up with anything better. Given that it is rather long, how would you feel about InfinitePF, similar to InfinitePEPS? (Along with PFTensor etc)

I don't really like the 'PF' acronym, for the simple reason that no one really uses it so there's no way of knowing what it means without actually reading the docstring. This is quite different for 'PEPS' I think. I'd be happy to use the shorter names as internal aliases, but I wouldn't like to export InfinitePF like that. Makes sense?

[pbrehmer]

Thanks a lot for the effort! To me the PR mostly looks quite nice already.

On the topic of naming things: I am a bit confused by calling the AbstractTensorMap{S,2,2} a partition function at all. Wouldn't a more generic name be MPOTensor (and InfiniteMPO)? There could probably be applications where you would want to contract or optimize a generic MPO which doesn't have to represent a partition function physically, so calling it PartitionFunctionTensor might be too specific. What do you think?

Also: let me know if I can contribute or finish up something! I anyway wanted to add 2D partition function support at some point so I'm happy to help.

[leburgel]

On the topic of naming things: I am a bit confused by calling the AbstractTensorMap{S,2,2} a partition function at all. Wouldn't a more generic name be MPOTensor (and InfiniteMPO)?

I agree using 'partition function' might imply some additional physical properties rather than just a network structure, which is really only what we mean here. The reason we didn't just reuse MPSKit.MPOTensor and MPSKit.InfiniteMPO is that these are really one-dimensional structures with some additional assumptions as well.

I would personally be okay with using the term 'partition function' for things with no stat mech origin, or at least more than I would be with using 'MPO' for 2D network structures. We could also go for SquareTensor and InfiniteGrid if we want to be more generic?

[pbrehmer]

Well, in principle partition functions could also be a 1D or 3D thing, right? So I feel a bit icky about using PartitionFunctionTensor in that way.

I do think it could make sense to use the MPSKit.MPOTensor alias because that does not yet imply a 1D structure. For InfiniteMPO I agree, we would need a different name (also to avoid a clash with MPSKit). What do you think about InfiniteMPOGrid or InfiniteMPONetwork?

[lkdvos]

To be fair, matrix product operator really does not make too much sense either in the 2D grid case...
The only other thing I can come up with is Projected Entangled Pair Function (PEPF), which might make some sense.
Ultimately, given that there is no term for this anywhere, we should just pick something that's practical for us though, PEPSKit.InfinitePartitionFunction does have the undertone of being 2D given that it is in PEPSKit. If you really want, we could also just add a type parameter N to give the dimensionality, which is reflected in the array type

[pbrehmer]

The InfiniteSquareNetwork interface already reduces the code duplication among InfinitePEPS, InfinitePEPO and InfinitePF quite a bit. I'm not yet sure, however, how to best incorporate the abstract type in sparse_environments.jl and ctmrg_contractions.jl. My plan would be to reduce the redundancies by having methods which could dispatch on either InfinitePF, Tuple{InfinitePEPS,InfinitePEPS} or eventually also Tuple{InfinitePEPS,InfinitePEPO,InfinitePEPS}. The goal is that only the methods that really contain the @tensor contractions specialize on the underlying tensor type - everything else (i.e. the methods that select the right coordinates) should be generic enough to only require one method per contraction, if that makes sense...

I'll finish up the rest tomorrow :-)

[pbrehmer]

There is a design decision to be made if we want to generalize our CTMRG code to also support different ket and bras, as well as PEPOs, and I'm not sure what the best way forward is. Right now, all CTMRG contractions support different ket and bra but there is no way for the user to pass this onto the CTMRG algorithm. And in principle our code is generic enough to take any state (which in principle should be a InfiniteSquareNetwork) which is then passed onto EnlargedCorner and the CTMRG contractions.

The problem is that there is no InfinitePEPSKetBra type which contains both the ket and the bra InfinitePEPS which the user can pass to leading_boundary. Similarly, if we want to support PEPOs, we would need a InfiniteKetPEPOBra type which contains the ket, bra as InfinitePEPSs and the InfinitePEPO. Should we implement such types (names are up for discussion :P) and subtype as InfiniteSquareNetwork or is that redundant somehow? Of course this is not strictly necessary to support 2D partition functions (the code runs as is) but I thought it would pay off to invest some time into generalizing this interface while we're already at it. I hope that makes sense!

[leburgel]

There is a design decision to be made if we want to generalize our CTMRG code to also support different ket and bras, as well as PEPOs, and I'm not sure what the best way forward is. Right now, all CTMRG contractions support different ket and bra but there is no way for the user to pass this onto the CTMRG algorithm. And in principle our code is generic enough to take any state (which in principle should be a InfiniteSquareNetwork) which is then passed onto EnlargedCorner and the CTMRG contractions.

I think it makes sense to have all the lowest level methods for things like leading_boundary and EnlargedCorner to work on a 'fully specified' network to be contracted. By this I mean either a partition function, a real pair of ket and bra PEPS, and so on. Passing just an InfiniteSquareNetwork to contains extra assumptions on the other layer in the network. This is fine, but it seems natural to make this a layer around the core routines, which always take fully specified disctinct layers. That also solver the issue of not being able to pass distinct ket and bra PEPS, since you could then do this by just skipping this first layer and directly calling the core routines.

The problem is that there is no InfinitePEPSKetBra type which contains both the ket and the bra InfinitePEPS which the user can pass to leading_boundary. Similarly, if we want to support PEPOs, we would need a InfiniteKetPEPOBra type which contains the ket, bra as InfinitePEPSs and the InfinitePEPO. Should we implement such types (names are up for discussion :P) and subtype as InfiniteSquareNetwork or is that redundant somehow?

I don't know if having these as actual types is necessary, it might be but it's not clear to me. Just a union type of tuples for all different possible sandwiches would do the job to start. I certainly don't think it makes sense to make these potential full sandwich types subtypes of InfiniteSquareNetwork, since then the hierarchy would become a bit circular I think. At least with the current meaning of InfiniteSquareNetwork these new types don't really naturally fall under it I think.

Of course this is not strictly necessary to support 2D partition functions (the code runs as is) but I thought it would pay off to invest some time into generalizing this interface while we're already at it. I hope that makes sense!

I'm a bit biased of course since I want to start using this as soon as I can, but purely in terms of compartmentalizing features I think it would make sense to just restrict this PR to working partition functions. Just to keep additions from growing bigger and bigger I would break out additional discussions as issues or follow up PRs. Of course I'm not in that much of a rush, but it seems we're moving away from the initial purpose a bit already.

[lkdvos]

I'm certainly not against that. I know we previously have done so, and called it PEPSSandwich, and PEPOSandwich. I do quite like those names 😉

[edit]
After reading @leburgel answer, I guess we really are straying away a bit too far from the PRs purpose, and are definitely cramming in way too many changes at once.

[pbrehmer]

I don't know if having these as actual types is necessary, it might be but it's not clear to me. Just a union type of tuples for all different possible sandwiches would do the job to start. I certainly don't think it makes sense to make these potential full sandwich types subtypes of InfiniteSquareNetwork, since then the hierarchy would become a bit circular I think. At least with the current meaning of InfiniteSquareNetwork these new types don't really naturally fall under it I think.

I have also considered this but the problem with tuples as sandwiches is that you need to be able to index the sandwich. But I do agree that it doesn't make sense to subtype the sandwiches as InfiniteSquareNetworks as well since it would create a circular hierarchy. Maybe it's fine to just add the corresponding indexing methods for tuples, e.g. Tuple{InfinitePEPS,InfinitePEPS}, but I'm not sure it will be pretty.

Anyways, I agree that we're putting too many additions into this PR and should probably generalize the CTMRG interface in another PR - I will open up an issue! I will finish up the contraction business and then the PR is good to go for me.

[pbrehmer]

Alright, I'd be fine with merging this as is :-) (once the lights turn green)

[leburgel]

I added an expectation_value method which takes a pair with a single index and a single rank-4 tensor map. While the updated syntax in the test case worked due to the added AbstractTensorMap{S,M,M} support, it seemed confusing to me to mix a tuple of indices with a single rank-4 tensor map. It seemed natural to have the simplest and most common case to work as intuitively as possible.

Aside from that, while the generalization of the local tensor contraction is nice, it's not covered by the tests and to be honest I can't immediately think of a relevant test withing the current test model. It's also not really clear to me what the inds are supposed to be in case of a single higher-rank tensor map, nor are there checks on the tensor map rank versus the grid size like there is for the matrix case.

I at least added a docstring to indicate what the current possibilities are. I guess since there's no concrete use case for the more general methods right now this is fine for the time being.

[pbrehmer]

Aside from that, while the generalization of the local tensor contraction is nice, it's not covered by the tests and to be honest I can't immediately think of a relevant test withing the current test model.

That's true and I'd like to add tests. Locally, I tested that at least the code runs and captures all the cases I could think of. To test for correctness, it would be nice to at least compare to an Ising correlation function perhaps? I'm not sure if there are exact expressions/integrals we could reference against.

It's also not really clear to me what the inds are supposed to be in case of a single higher-rank tensor map, nor are there checks on the tensor map rank versus the grid size like there is for the matrix case.

You do need the inds since a higher-rank tensor map does not contain any information on how it is contracted with the surrounding environment. To make an example, a AbstractTensorMap{S,4,4} could be a 2x2 plaquette or it could be a 3x1 string. The first case would be specified using inds=((1, 1), (2, 1), (1, 2), (2, 2)) whereas the second case would correspond to inds=((1, 1), (2, 1), (3, 1)). That also shows why a general check for the tensor map rank is not as easy as for the matrix case since the number of supplied indices could technically vary given a fixed tensor rank.

[leburgel]

Locally, I tested that at least the code runs and captures all the cases I could think of. To test for correctness, it would be nice to at least compare to an Ising correlation function perhaps? I'm not sure if there are exact expressions/integrals we could reference against.

To just check the code I think we could use a 2x2 patch with one magnetization/energy tensor and the other three just regular partition function tensors, both in matrix form and blocked to cover all cases.

You do need the inds since a higher-rank tensor map does not contain any information on how it is contracted with the surrounding environment. To make an example, a AbstractTensorMap{S,4,4} could be a 2x2 plaquette or it could be a 3x1 string. The first case would be specified using inds=((1, 1), (2, 1), (1, 2), (2, 2)) whereas the second case would correspond to inds=((1, 1), (2, 1), (3, 1)). That also shows why a general check for the tensor map rank is not as easy as for the matrix case since the number of supplied indices could technically vary given a fixed tensor rank.

I see. Then the docstring I added should cover the proper usage I think, and maybe we could add a consistency check on the rank after determining the patch size from the given inds? At that point it's clear what the compatible tensor rank should be, right? At least this would be better than to let it throw an error for a badly specified contraction.

[pbrehmer]

To just check the code I think we could use a 2x2 patch with one magnetization/energy tensor and the other three just regular partition function tensors, both in matrix form and blocked to cover all cases.

Then the docstring I added should cover the proper usage I think, and maybe we could add a consistency check on the rank after determining the patch size from the given inds? At that point it's clear what the compatible tensor rank should be, right? At least this would be better than to let it throw an error for a badly specified contraction.

Good ideas! I will cook this up now.

[lkdvos]

I left a few remarks in the code, but honestly this is getting absolutely intractable to make sense of what is going on... (+1833,-519) lines of code, split over multiple topics is really hard to follow...
For future PRs, I would suggest we try and keep features and changes separated into different PRs, and before making large organizational overhauls, make sure that everyone is on board with the changes?

Let's try and simplify this so we can merge this ASAP, and then re-evaluate what needs to change. I would suggest the following:

1. This PR is mostly about having CTMRG for partition functions, so I would really refrain from going overboard with the features such as elaborate expectation_value variants, which we can definitely add in a separate PR. Can we maybe fix the argument order to have the interface in line with MPSKit, and simply have a working implementation of just inserting an operator into the partition function? We can always refactor later.
2. Given that this is now intimately tied to the InfiniteSquareNetwork rewrite, I'm okay with merging that as is, but can we re-discuss how we want to organize that here: #116, and then decide on an implementation afterwards?
3. Can we open a separate issue/PR for a round of docstring fixes where we first decide on a style, and then work together to implement that?

[lkdvos]

This sounds to me like something that is not too well defined. It's a bit odd to me to now all of a sudden use inds to signify rectangular regions instead of actual positions, and since we don't have any conventions for how this plaquette tensor would then look like, why are we defining this? Can we simply do this later?

[pbrehmer]

Yes, there is an ambiguity in working with higher-rank tensor maps and supplying the inds. The annoying thing in the current implementation is that only really the indices that specify the corners of the rectangular patch matter, so I agree that it's not a nice way of implementing multi-site operators.

Maybe to restore the original way of using the inds we could also do the following: Instead of supplying a higher-rank tensor map or matrix of tensor maps, we could just supply a tuple of Pair{CartesianIndex(2),TensorMap{S,2,2}} (similar to how a LocalOperator works) as well as the partition function to the contract_local_tensor function. Then the supplied tensors would just replace the partition function on the supplied indices. (If no tensors are supplied, then the function just computes the norm as in contract_local_norm.)

I don't want to further blow up the PR so I'm also happy to go for a simpler alternative. But I think the main contractions are already implemented to generalize the partition function operator evaluation to it should just be a matter of reorganizing that.

[leburgel]

Yes, there is an ambiguity in working with higher-rank tensor maps and supplying the inds. The annoying thing in the current implementation is that only really the indices that specify the corners of the rectangular patch matter, so I agree that it's not a nice way of implementing multi-site operators.

For now, I would just remove the higher-rank tensor case. It's a bit random to come up with an interface for this kind of thing if we don't even have a use case in mind. Since this is making things more difficult here, I suggest we drop it altogether and think about reintroducing it once we have a specific need for it.

Maybe to restore the original way of using the inds we could also do the following: Instead of supplying a higher-rank tensor map or matrix of tensor maps, we could just supply a tuple of Pair{CartesianIndex(2),TensorMap{S,2,2}} (similar to how a LocalOperator works) as well as the partition function to the contract_local_tensor function. Then the supplied tensors would just replace the partition function on the supplied indices. (If no tensors are supplied, then the function just computes the norm as in contract_local_norm.)

This seems reasonable, and also allows to solve the above issues. I don't really see the need to have a version of contract_local_tensor where no tensors are supplied though. It would be more transparent to just do the explicit call where you pass the partition function tensors (I assume you're referring to the denominator in expectation_value?). A rename to contract_local_tensors would then be appropriate of course...

[pbrehmer]

I would definitely agree that we should scrap the AbstractTensorMap{S,M,M} case since that is not nice to work with and to implement.

I could still see that the rest is salvageable :-) I don't find the case of correlation functions too random - so if someone wants to evaluate e.g. nearest-neighbor correlators that is not super exotic, right? If you want I could still try to clean up/simplify the general contraction mechanism using Vector{Pair{CartesianIndex{2},AbstractTensorMap{S,2,2}}} so that it is more in line with the LocalOperator way of computing the expectation value. But if we want to scrap that for now, I'd also be fine with that - so feel free to remove the relevant code bits.

[lkdvos]

In any case, this is completely decoupled from this PR, so we can separate that out to get reorganized :)

[pbrehmer]

Yes let's do that. Then let's remove the multi-site operator functionality from this PR (reverting back to @leburgel original 1-site operator code) and then I will open up another PR where we can discuss this further :-) In there I can push my current idea of what would be the most organized approach with the least ambiguities.

[lkdvos]

Before you go all-in with implementing, we might want to discuss first, unless you are okay with possibly having to make large changes if we find a better approach?

[leburgel]

Yes let's do that. Then let's remove the multi-site operator functionality from this PR (reverting back to @leburgel original 1-site operator code) and then I will open up another PR where we can discuss this further :-) In there I can push my current idea of what would be the most organized approach with the least ambiguities.

I reverted (not a clean one, but everything is still in the commit history here anyway) back to the simple one-site case and matched the expectation_value syntax to the other methods. Happy to discuss extensions in follow-ups. I'll try to remember what kinds of stat mech quantities we might want to be able to compute, it's been a while :)

[pbrehmer]

Before you go all-in with implementing, we might want to discuss first, unless you are okay with possibly having to make large changes if we find a better approach?

Sure, I definitely had discussions in mind. What I was planning could maybe serve as a good starting point for discussion/further refinements.

[lkdvos]

I think we can merge this like this if the tests all pass.