Generalize InfinitePartitionFunction observable evaluation

[pbrehmer]

This PR will add generic contractions which allow us to use multi-site operators in expectation_value(pf::InfinitePartitionFunction, ...). Here, we will try to reuse the local operator contraction mechanism which automatically constructs the correct enlarged environment and inserts states/operators based on the supplied inds.

[codecov]

Codecov Report

Attention: Patch coverage is 94.82759% with 3 lines in your changes missing coverage. Please review.

Files with missing lines	Patch %	Lines
src/algorithms/contractions/localoperator.jl	94.44%	3 Missing ⚠️

Files with missing lines	Coverage Δ
src/algorithms/toolbox.jl	97.50% <100.00%> (+0.06%)	⬆️
src/algorithms/contractions/localoperator.jl	94.77% <94.44%> (-0.43%)	⬇️

[pbrehmer]

I have added an initial sketch of what I was talking about yesterday in #111. Operators that act on N sites are supplied via a NTuple{N,Pair{CartesianIndex(2),AbstractTensorMap{S,2,2}}}, so e.g. ops = ((1, 1) => op1, (2, 1) => op2) for a nearest-neighbor operator.

While I do think that this is functional, I'm not sure about a few things. For example, to compute the denominator in expectation_value, I reuse contract_local_tensors by supplying nothing instead of the operator - not sure if I like that.

Overall, I'd be okay with having an interface like this since you can readily evaluate correlation functions or arbitrary operator strings and it is kind of similar to what we have with LocalOperator with the difference that we only specify on-site tensors and can't have higher-rank tensors. But of course I'm not fixed to this idea and I'm happy to discuss and discard this approach if necessary.

[lkdvos]

The thing I don't like about this approach is that it is fully biased towards non-symmetric tensors, and does not generalize. For example, even in the Ising model, inserting a single flip in a symmetric partition function does not work with the Z2 symmetry, and really what you need is an auxiliary leg to carry that charge, which connects the other flip, in order to obtain a non-zero result. In that sense, this is not really a clean generalization of the quantum case, where you can insert pairs of spin flips at 2 locations because they are represented by a single operator. Additionally, the syntax in terms of ntuples doesn't provide enough information, because even if we would give each tensor an additional leg, you now also need to specify how to pair them up in the case where N>2.
In general, I would say that this is the kind of code that is a bit against the design of TensorKit, given that it is not agnostic about the symmetry. If at all possible, I would really prefer to have an API that does not need to know about the symmetry first, before considering the implementation.

Given these constraints, I would maybe suggest that we don't try to push two-point correlators into the same formalism, given that there is not really that much shared code anyways. In MPSKit.jl, there is a function

"""
    correlator(ψ, O1, O2, i, j)
    correlator(ψ, O12, i, j)

Compute the 2-point correlator <ψ|O1[i]O2[j]|ψ> for inserting `O1` at `i` and `O2` at `j`.
Also accepts ranges for `j`.
"""

while I don't think the second syntax makes sense here, I would be more tempted to add the first. This is a more intuitive syntax for what we're trying to compute anyways, and does not need to include the generic case of N>2, so we could have O1::AbstractTensorMap{T,S,2,3} being "excited" partition function tensors, which can now uniquely be matched up.

[pbrehmer]

Indeed, I forgot about symmetric tensors. Then let's forget about the N>2 case and try to implement the correlator approach as in MPSKit. I would have two questions:

• Is it necessary that O1 and O2 are adjacent, i.e. nearest-neighbor tensors? Or can the auxiliary leg also be connected across multiple sites?
• Should we still allow for AbstractTensorMap{T,S,2,2} as an input (since that is the usual case of non-symmetric tensors) and then just convert to an AbstractTensorMap{T,S,2,3} (or AbstractTensorMap{T,S,3,2}) with a dummy leg?

[lkdvos]

It's definitely not necessary for both to be adjacent, and I have no problem also supporting the case without an additional leg. The bigger question is whether the implementation via a full contraction is the correct one. Keep in mind that this code is written for evaluating local expectation values, ie nearest-neighbour, next-nearest neighbour etc, and really does not scale very well. On the one hand, if the correlator would be in a vertical or horizontal line, you really would not want to hand off this entire contraction to @tensor with @autoopt, because the cost of that scales exponentially with the number of tensors, while we actually now how to contract that. For diagonal correlators, there are probably smarter approaches using channel methods etc. It's not because we can write down the contraction for such a correlation function that we can actually contract it 😉

[pbrehmer]

Well, you have a point. Probably the best way to go forward is to wait until someone actually needs to evaluate correlators and then implement it (well who could have thought :P). Since I'm short on time and really need to get some other stuff done, I'll probably let the PR sleep for a while or someone else implements this.