[Yue-Zhengyuan]

Exact contraction becomes costly for observables spanning large lattice patches.
This PR adds boundary-MPS contraction for single-layer PEPO expectation values through an algorithm-dispatched expectation_value API.

expectation_value(ρ, O, env) # ExactContract() remains the default

alg = BoundaryMPSContract(; trunc = truncrank(χ), maxiter = 1, direction = :auto)
expectation_value(ρ, O, env, alg)

• Supports dense, MPO, and tensor-product terms through the existing LocalOperator interface.
• Routes observables through finite patches and normalizes each term independently.
• Uses zipup truncation (Zip-up compression for FiniteMPO-FiniteMPS product MPSKit.jl#470) with optional DMRG refinement and automatic or explicit row/column sweeps.
• Shares term evaluation through algorithm dispatch, with exact single-site contractions.
• By default, the truncation rank follows the largest CTMRG boundary dimension.

[codecov]

Codecov Report

❌ Patch coverage is 96.56085% with 13 lines in your changes missing coverage. Please review.

Files with missing lines	Patch %	Lines
...rc/algorithms/contractions/mpo_path/pepo_1layer.jl	92.79%	8 Missing ⚠️
src/algorithms/contractions/mpo_path/routing.jl	98.23%	2 Missing ⚠️
src/algorithms/expectation_value/boundary_mps.jl	88.88%	2 Missing ⚠️
src/algorithms/contractions/patch/pepo_1layer.jl	98.70%	1 Missing ⚠️

Files with missing lines	Coverage Δ
src/PEPSKit.jl	100.00% <ø> (ø)
src/algorithms/contractions/patch/tools.jl	100.00% <100.00%> (ø)
src/algorithms/contractions/transfer.jl	100.00% <100.00%> (+4.34%)	⬆️
src/algorithms/expectation_value/algorithms.jl	100.00% <100.00%> (ø)
.../algorithms/expectation_value/expectation_value.jl	100.00% <100.00%> (ø)
src/utility/indexing.jl	92.50% <ø> (ø)
src/algorithms/contractions/patch/pepo_1layer.jl	98.70% <98.70%> (ø)
src/algorithms/contractions/mpo_path/routing.jl	98.23% <98.23%> (ø)
src/algorithms/expectation_value/boundary_mps.jl	88.88% <88.88%> (ø)
...rc/algorithms/contractions/mpo_path/pepo_1layer.jl	92.79% <92.79%> (ø)

... and 4 files 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]

@lkdvos @leburgel Does the MPSKit support for fermions (via planar contraction) still work if the MPS/MPO arrows do not follow the standard directions?

[leburgel]

@lkdvos @leburgel Does the MPSKit support for fermions (via planar contraction) still work if the MPS/MPO arrows do not follow the standard directions?

Yes, I think it should just work regardless of which spaces exactly are dual.

[Yue-Zhengyuan]

Fermion support is now here, but the way to achieve it is not ideal. I first put the virtual arrow directions in the iPEPO in the standard direction. But the on the right boundary, the "physical space" is dual, as required by the CTMRGEnv constructor (see C₂ below):

    [1; 2]      [1 2; 3]        [1; 2]
    C₁-←-2      1-←-E₁-←-3      1-←-C₂
    ↓               ↓               ↑
    1               2               2

Thus I still need to use many planar operations (some of which I still don't fully understand) when constructing the left and right ends of each row, and the south boundary.

(Edit: now I also flip the east envspace to standardize all arrows. Things then make more sense.)

Otherwise, if I dispatch away from planar operations, it seems that I need to overload lots of things in MPSKit (not only in ZipUp but also DMRG for further refinement), which may not be worthwhile. (But anyway we still need to do it when generalizing to approximate contraction of multi-layer networks?)

[leburgel]

I was going to take some time to go over this, but it seems this PR is rather large and I completely misjudged the time I allocated. I'll try to go over the details as soon as I can, but maybe I can start with a few general remarks already before I submit a full review:

---

As far as I can tell, the newly introduced MPOObservable could in spirit serve as a term in a LocalOperator, where instead of associating a dense TensorMap to a group of sites we associate a finite MPO (encoded here as a vector of MPO tensors). The only difference is the path field that is additionally included in the MPOObservable struct. This leads me to the question: what is the benefit of being able to manually specify a path, over just determining this automatically? Also, what exactly is the difficulty in this automatic path finding (as indicated in the note above)? If it's not more expensive to always use a default path, then we can directly use MPOs as LocalOperator terms, whose expectation value we could then evaluate exactly or approximately. This would be very nice to have, hence my question.

---

For the finite window contraction, would it be possible to directly pass an algorithm for this? I can imagine wanting to just use DMRG2 (or even just DMRG) to do this, rather than starting from ZipUp and then optionally refining. Probably the way it is now is fine for practical purposes, at least the ones you had in mind, but it feels like it would be useful to have direct control over the algorithm.

---

Regarding fusing the MPO path, instead of fusing everything into the virtual spaces, we could (at some point) try to keep the virtual MPO legs around explicitly as auxiliary legs for MPS tensors when they are encountered. The problem is then to contract two "matching" auxiliary legs corresponding to the same virtual MPO bond when they encounter each other in two tensors being contracted. This is a lot more tedious to do in our framework, since we don't have "labeled" tensor indices that recognize each other automatically. Here we can just pass around explicit labels for all of the auxilary legs that come from a virtual MPO bond, and keep track of them. Not sure how feasible this is, and definitely not sure if this is something we want to do here, but just wanted to bring it up already.

---

Regarding fermion support, I'm not sure if I'm entirely following the discussion. I think the MPSKit routines should just be able to handle physical spaces with a "non-standard" duality. Using planar operations where possible should ensure this just works, in the way that is consistent with MPSKit conventions. Therefore, I don't really see the need to consider dispatching away from planar operations. The only reason for this as far as I'm concerned would be to ensure differentiability, which is not a goal at all for now I think. So even flipping this east env space shouldn't be necessary, I think.

The only reason to not use planar operations should be if we simply can't, e.g. for multi-layer networks. There, we have to overload the core contractions used in MPSKit anyway, but we additionally have to put additional twists to ensure the result is what is expected by MPSKit conventions. For the MPSKit contractions already overloaded in PEPSKit, this should have been done in a way that they can also handle arbitrary arrow directions.

---

For multiplying in the bottom boundary, I would really like to avoid (doubly) conjugating edge tensors, this always leads to issues sooner or later in my experience. The question then becomes if we need to twist anything on the resulting top MPS (after all of the MPO applications) to ensure we can just contract in south edges in the usual way and actually get what we want. I'm not entire sure about that.

[Yue-Zhengyuan]

In principle I should split this PR in two: first focusing only on expectation_value_approx, with a follow up with correlator_approx. But these two share much infrastructure, so I put them together...

• Path routing for MPO terms

I agree that needing an additional path to specify the MPO is not ideal. The difficulty is in writing a clever algorithm that creates a path that go through the MPO sites in order, while also being non-self-intersecting. You can consider the sites (1, 1) - (6, 1) - (2, 0) - (3, 2) - (4, 0) - (5, 2) (which is unlikely to appear in practice, though). It is not obvious to me which should be the "canonical" path joining them without self-intersection.

• Algorithm struct for window contraction

It is already there, called WindowApprox. The problem with skipping ZipUp is how you initialize a ψ0 for DMRG/DMRG2 without randomly guessing.

Currently this struct is not exported. Instead, we provide kwargs trunc and maxiter in the API to build the struct internally:

WindowApprox(Zipup(; trunc), _approx_dmrg(maxiter))

• Fusing MPO virtual space with the PEPO

I want to postpone this. Currently ZipUp (and DMRG?) only support(s) MPS with one physical leg per site.

• Fermion support

Changing arrow directions in MPSKit is like planar contraction with flippers without twists. But this is not the case in PEPSKit, where we should use flip that preserves @tensor calls. First flipping all arrows to the standard direction appears to automatically put the required twists back in the network; the change of arrows is more like a by-product that is unimportant to MPSKit. I suggest we keep this pre-processing for now, and only revisit it at the end.

• Conjugation of south boundary

In the cache WindowRowCache used by correlator_approx, we not only need the south CTM boundary, but also need the contraction of the south boundary with some rows above it. So avoiding conjugating south boundary twice will require more than a custom version of dot. We will also need to write ψ * O, where the MPO O acts on the bra-state ψ. This is a bit redundant to me.

[leburgel]

Sorry it took me so long to get through this, I couldn't seem to find the time to sit down and get through the whole thing in one go.

I have a first round of comments, most of which are related to the general remarks I had a while back. An overview:

• Path routing for MPO terms: I think we can make this problem much easier by giving up on ordering the sites according to the matrix-like linear ordering we use for dense local operators. There the choice doesn't matter much so imposing a canonical ordering makes sense. But for MPOs, choosing the site ordering in a way that makes the routing easy makes much more sense to me.

• Algorithm struct for window contraction: I left a more in-detail comment below, but I think that if we use DMRG with an expansion we can merge things quite cleanly actually.

• Fusing MPO virtual space with the PEPO: definitely fair, fusing is the quickest way to actually get something up and running.

• Fermion support: as always I am still confused by this. If it's verified to work properly, there's not many objections I can raise, so that would be fine for me. We can always revisit how we do things later, but I don't have any conceptual arguments right now.

• Conjugation of the south boundary: again there's a more detailed comment below, but I think we can avoid conjugations by rotating instead. This comes with its own issues, but I think it's conceptually cleaner than taking adjoints. To implement ψ * O, we would just rotate O, compute the normal thing, and take care to reverse ψ back and forth in the appropriate way. This is the same way I would do infinite boundary MPS contractions, so I think that could be a viable alternative.

[leburgel]

I would personally implement this through rotation without conjugation, which reverses the left-right ordering of the south edge tensors and collects them into a finite MPS. This respects the natural virtual orientation of the south edge tensors, and avoids having to take the "adjoint" of individual MPS tensors.

[Yue-Zhengyuan]

I think we still need to handle the physical legs if we use the rotation approach? South boundary will have dual physical spaces.

[leburgel]

I somehow thought dual physical spaces weren't really an issue for MPSKit, I though we only made hard assumptions on the virtual space duality. I might be wrong or be misremembering, I just thought I tried this out at some point a long time ago.

[leburgel]

The main reason why I originally asked if we could get rid of this path field is that this would then just become a specific term in a LocalOperator. Conceptually, I quite like just being able to call expectation_value_approx on any LocalOperator, where we then evaluate every term individually in the appropriate way (using some kind of local_expectation_value_aproximate hook).

For dense tensors, we could just convert them to an MPO when calling the approximate version. The path routing could then happen inside this local expectation value computation, concerning only the relevant path.

The non-intersecting path routing problem is indeed hard in general, if we're handed an MPO acting on specific sites. But usually, we would want to evaluate some operator which we can represent as a dense tensor map (I assume this is the most common case?), which we then decompose into an MPO and evaluate approximately. If we choose the site order in this decomposition in a way that takes into account the routing that we'll have to do, this makes the problem much easier. So I think we could get rid of the routing difficulties by simply allowing sites that are not ordered according to the lattice ordering. This anyway makes a lot of sense for MPO observables I think, where intuitively we would order the sites to end up with the virtual string of smallest length, which will generally be quite easy to route.

Does this make sense, and sound like a reasonable thing to try?

In the end I would like to not do any path routing and virtual fusing at all, and just pass around the MPO indices appropriately automatically. This is not at all necessary now, but it would be good to already get ahead of the routing issue.

[Yue-Zhengyuan]

For expectation_value, we may add an additional algorithm struct like ExactContract for the existing exact contraction, and dispatch to boundary MPS approximation using the contraction algorithm:

expectation_value(state, O, env, [alg = ExactContract()])
expectation_value(state, O, env, alg = WindowApprox(zipup, dmrg))

So extra names like expectation_value_approx will be unnecessary.

correlator_approx may still be kept, since it accepts a list of 2-site bonds instead a fixed first site + a list of second sites. (But I may also make the transition?)

I think even in common cases, the MPO is not always first constructed as a dense map? (e.g. some are just a TensorProductTerm) And it is not easy either to come up with a algorithm that reorder the sites acted on by the operator so that the string length is minimized.

[leburgel]

If we let go of the operator site ordering being the same as for dense tensor local terms, I think we could simplify things quite a bit?

[Yue-Zhengyuan]

Major changes made today (mainly on correlator_approx):

• Changed correlator_approx to accept one fixed site i and second sites js, matching correlator. Direction is selected automatically from the target positions; square windows prefer row-by-row contraction.
• Introduced adaptive window depth with fixed width: each contraction ends at its target row or column.
• Removed WindowRowCache, the reverse sweep that needs _adjoint_mpo, and operator swapping.

[Yue-Zhengyuan]

Now the only big issue is designing the MPO terms. Ideally this should also include PBC-MPOs that form a closed loop. I'll take a closer look at #425 and see if I can get some inspirations.

[Yue-Zhengyuan]

I'll move correlator_approx elsewhere to make this PR focused on expectation_value.